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Abstract 

In this paper we present a new approach to modeling finite set domain constraint prob- 
lems using Reduced Ordered Binary Decision Diagrams (ROBDDs). We show that it is 
possible to construct an efficient set domain propagator which compactly represents many 
set domains and set constraints using ROBDDs. We demonstrate that the ROBDD-based 
approach provides unprecedented flexibility in modeling constraint satisfaction problems, 
leading to performance improvements. We also show that the ROBDD-based modeling 
approach can be extended to the modeling of integer and multiset constraint problems in 
a straightforward manner. Since domain propagation is not always practical, we also show 
how to incorporate less strict consistency notions into the ROBDD framework, such as set 
bounds, cardinality bounds and lexicographic bounds consistency. Finally, we present ex- 
perimental results that demonstrate the ROBDD-based solver performs better than various 
more conventional constraint solvers on several standard set constraint problems. 



1. Introduction 

It is often natural to express a constraint satisfaction problem (CSP) using finite domain 
variables and relations between those variables, where the values for each variable are drawn 
from a finite universe of possible values. One of the most common methods of solving finite 
domain CSPs is through maintaining and updating a domain for each variable using the 
combination of a backtracking search and an incomplete local propagation algorithm. The 
local propagation algorithm attempts to enforce consistency on the values in the variable 
domains by removing values that cannot form part of a solution to the system of constraints. 
Various levels of consistency can be defined, each with associated costs and benefits. Most 
consistency algorithms are incomplete — that is, they are incapable of solving the problem 
by themselves, so they must be combined with a backtracking search procedure to produce 
a complete constraint solver. 

When attempting to apply this general scheme to the task of solving constraint satis- 
faction problems over finite set variables we quickly run into practical problems. Since the 
universe of possible values for a set variable is usually very large, the naive representation of 
the domain of a set variable as a set of sets is too unwieldy to solve realistic problems. For 
example, if a set variable can take on the value of any subset of the set {1, . . . , N}, then its 
domain contains 2^ elements, which quickly becomes infeasible to represent as N increases 
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in magnitude. Accordingly, most set constraint solvers to date have used an approximation 
to the true domain of a set variable in order to avoid this combinatorial explosion. 

The most common approximation to the true domain of set variable v has been to 
maintain an upper bound U and a lower bound L of the domain under the subset partial 
ordering relation and to perform set bounds propagation on these bounds. That is, L 
contains elements that must be in the set v, and U is the complement of the set of elements 
that must not be in v. Conventionally, a fixed set of inference rules specific to each constraint 
are used to enforce consistency on these upper and lower bounds. This basic scheme was 
proposed by Puget (1992), and has been implemented by set solvers such as Conjunto 
(Gervet, 1997), the fd.sets and ic.sets libraries of ECLTS 6 (IC-PARC, 2003), Mozart 
(Miiller, 2001), and ILOG Solver (ILOG, 2004). 

Set bounds are a crude approximation to a set domain at best, and thus various re- 
finements to the bounds propagation scheme have been proposed so as to more effectively 
capture the nature of a set domain. Azevedo (2002) demonstrated that maintaining and 
performing inferences upon the upper and lower bounds on the cardinality of a set domain 
leads to a significant performance improvement on a variety of problems. While earlier 
set solvers such as Conjunto also maintained cardinality bounds, only partial usage was 
made of this information. More recently, Sadler and Gervet (2004) showed that incorporat- 
ing upper and lower bounds under a lexicographic ordering leads to significantly stronger 
propagation, albeit at the cost of a marked increase in propagation time, leading to only a 
marginal performance improvement overall. While both of these approaches provide more 
effective propagation than the simple set bounds scheme, they do not approach the effec- 
tiveness of a true set domain propagator, which ensures that every value in the domain 
of a set variable can be extended to a complete assignment of every variable in any given 
constraint. 

Consequently, we would like to devise a representation for set domains and constraints 
on those domains that is tractable enough to permit domain propagation. The principal 
observation that permits the implementation of a set domain propagator is that a finite 
integer set v can be represented by its characteristic function \v'- 

Xv '■ Z — > {0, 1} where Xvif) = 1 iff i € u 

Accordingly we can use a set of Boolean variables V{ to represent the set v, which correspond 
to the propositions Vi <-> i G v. We can describe set domains and set constraints in terms 
of these Boolean variables. Interestingly, set bounds propagation as described above is 
equivalent to performing domain propagation in a naive way on this Boolean representation. 

In this paper we investigate this Boolean modeling approach for modeling finite domain 
constraints, and in particular set constraints. We show that it is possible to represent the 
domains of set variables using Reduced Ordered Binary Decision Diagrams (ROBDDs), a 
data structure for representing and manipulating Boolean formulae. Such representations 
are usually fairly compact, even if they correspond to very large domains. In addition, it 
is possible to represent the constraints themselves as ROBDDs, permitting us to produce 
efficient set domain constraint propagators solely using ROBDD operations. 

The ROBDD-based representation allows us to easily conjoin constraints and existen- 
tially quantify variables, thus permitting us to remove intermediate variables and merge 
constraints for stronger propagation. The construction of global constraints becomes an 
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almost trivial exercise, without the requirement to write laboriously new propagators for 
every new constraint that we would like to use. 

We also demonstrate that only minor changes are needed to the operation of the set 
domain propagator in order to implement other, less strict notions of consistency. In par- 
ticular, we show how to construct set bounds, cardinality bounds and lexicographic bounds 
propagators, and how to utilize these notions to produce a split domain solver that combines 
bounds and domain reasoning. 

A key theme of this paper is the flexibility of the ROBDD-based modeling approach. 
We are not limited only to modeling set variables — for example ROBDDs can be used 
to model integer variables and integer constraints. While the ROBDD-based approach 
to modeling finite domain integer variables is in general not as efficient as many existing 
finite domain integer solvers, the ability to model integer constraints using ROBDDs is an 
essential building block which allows us to construct set constraints such as cardinality and 
weighted sum, as well as allowing us to represent multisets and multiset constraints. 

Finally, we present experiments using a variety of standard set problems to demonstrate 
the advantages of our modeling approach. 

Many of the ideas from this paper have been previously published in two previous works 
(Lagoon & Stuckey, 2004; Hawkins, Lagoon, & Stuckey, 2004). This paper contains a more 
complete exposition of those ideas, as well as important extensions of the work that has 
previously been presented. These include substantial improvements in the complexity of the 
propagation algorithm, as well as cardinality and lexicographic bounds solvers implemented 
using ROBDDs. In addition, we show how to model integer expressions, allowing us to 
implement a weighted-sum constraint, and to model multisets and multiset constraints. 
Using this, we present new experimental results for the Hamming Code and Balanced 
Academic Curriculum problems. Finally, we present results comparing the ROBDD-based 
solver against a solver with good cardinality reasoning (Mozart). 

The remainder of this paper is structured as follows. Section 2 contains essential concepts 
and definitions necessary when discussing finite domain solvers and ROBDDs. Section 3 
shows how to model set domains and set constraints using ROBDDs, and presents a basic 
outline of an ROBDD-based set constraint solver. Section 4 demonstrates how to improve 
the performance of the ROBDD-based set solver through the construction of global con- 
straints, the removal of intermediate variables, and through symmetry-breaking approaches. 
Section 5 demonstrates how to model integer and multiset expressions, as well as how to 
implement a weighted-sum constraint for set and multiset variables. Section 6 shows how 
to construct a more efficient domain propagator, as well as set bounds, set cardinality, and 
lexicographic bounds propagators. Finally, in Section 7 we present experimental results 
comparing the ROBDD-based solver with more conventional set solvers on a variety of 
standard problems. 

2. Preliminaries 

In this section we define the concepts and notation necessary when discussing propagation- 
based constraint solvers. These definitions are largely identical to those presented by Lagoon 
and Stuckey (2004) and others. For simplicity we shall present all of our definitions in the 
case of finite set variables; the extensions to multiset and integer variables are trivial. 
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2.1 Lattices, Domains, and Valuations 

Let C be the powerset lattice (P(U),c), where Vix) denotes the powerset of x and the 
universe U is a finite subset of Z. We say a subset M C £ is convex if for any a,c £ M 
the relation aGCc implies b G M for all b £ C. The interval [a, b] is the set M = {x G 
£ I a C x C 6}. Intervals are obviously convex. Given a subset K C £, we define the convex 
closure of if: 



The convex closure operation satisfies the properties of extension (x C com; (a;)), idempo- 
tence (conv(x) = conv(conv(x))), and monotonicity (if x C y then conv(x) C (y)) (Gervet, 



Example 2.1. The set X = {{1}, {1, 3}, {1,4}, {1, 3, 4}} is convex and equivalent to the in- 
terval [{1}, {1, 3, 4}]. Conversely, the set Y = {{1}, {1, 3}, {1, 3, 4}} is not convex. However, 
the convex closure of Y is precisely the interval X, i.e. conv(Y) = X. 

Let V denote the fixed finite collection of all set variables. Each variable has a domain, 
which is a finite collection of possible values from £ (which are themselves sets). More 
generally, we define a domain D to be a complete mapping from the set V to finite collections 
of finite sets of integers. When we speak of the domain of a variable v, we mean D(v). A 
domain D± is said to be stronger than a domain D2, written D\ C D2, if Di(v) C D2(v) 
for all v G V. Two domains are said to be equal, written D\ = D2, if D\(v) = D2(v) 
for all v G V. We call a domain D a range domain if D(v ) is an interval for all v G V. 
We extend the concept of convex closure to domains by defining ran(D) to be the unique 
(range) domain such that ran(D)(v) = conv(D(v)) for all dGV, 

A valuation is a function from V to £, which we write using the mapping notation 
{vi 1 ^ di, -U2 1— ► c?2, • • • , v n 1— > d n }, where € V and di £ C for all i 6 N. Clearly a valuation 
can be extended to constraints involving the variables in the obvious way. We define vars to 
be the function that returns the set of variables that are involved in a constraint, expression 
or valuation. In an abuse of notation, a valuation is said to be an element of a domain 
D, written 9 G D, if 0(v) € D(v) for all v G vars{6). 

We say a domain D is a singleton or valuation domain if |-D(f)| = 1 for all v G V. In 
this case D corresponds to a single valuation 6 where D{v) = {0(v)} for all d€V. 

2.2 Constraints, Propagators, and Propagation Solvers 

A constraint is a restriction placed upon the allowable values for a collection of variables. 
We are interested in the following primitive set constraints, where k is an integer, d is a 
ground set value, and n, v and w are set variables: k G t> (membership), k ^ v (non- 
membership), u = v (equality), u = d (constant equality), u Q v (non-strict subset), 
u = v U w (union), u = vPiw (intersection), u = v\w (set difference), u = v (complement), 
u ^ v (disequality), \u\ = k (cardinality), |u| > k (lower cardinality bound), and \u\ < k 
(upper cardinality bound) . Later we shall introduce non-primitive set constraints which are 
formed by composing primitive set constraints. 

We define the solutions of a constraint c to be the set of valuations that make that 
constraint true, i.e. solns(c) = {6 \ vars(9) = vars(c) and N 9(c)}. 
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Example 2.2. Suppose v and w are set variables, and D is a domain with D(v) = 
{{1}, {1,3}, {2, 3}}, and D(w) = {{2}, {1, 2}, {1, 3}}. Let c be the constraint v C w. 
Then the solutions to c in the domain D are the valuations {v i— > {l},iu ► {1,2}}, 
{f i ^ {l},w | — ► {1,3}}, and {v t— > {1,3}, it; i— > {1,3}}. 

With every constraint we associate a propagator f, which is a monotonic decreasing 
function from domains to domains, so if D\ C Z?2 then f(D\) C /(-D2) and f(D\) C £>i. A 
propagator / is said to be correct for a constraint c if: 



Correctness is not a strong restriction, since the identity propagator is correct for all con- 
straints c. We usually assume that all propagators are both correct and checking, that is, if 
D is a singleton domain formed through propagation then its associated valuation 9 makes 
9(c) true. 

We can use propagators to form the basis for a constraint solver. A propagation solver 
solv(F,D) takes a set of propagators F and a current domain D, and repeatedly applies 
propagators from F to the current domain until a fixpoint is reached. In other words 
solv(F,D) is the weakest domain D' C D which is a fixpoint (i.e. f(D') = D') for all 
/ G F. This fixpoint is unique (Apt, 1999). 

2.3 Local Consistency 

The notion of local consistency is of importance when considering the solution of constraint 
satisfaction problems. We can define various levels of consistency of different strengths and 
levels of difficulty to enforce. We describe two here. 

A domain D is said to be domain consistent for a constraint c if D is the strongest 
domain which contains all solutions 9 € D of c; in other words there does not exist a 
domain D' C D such that 9 G D and 9 G solns(c) implies 9 G D'. A set of propagators 
maintains domain consistency if solv(F,D) is domain consistent for all constraints c. 

Definition 1. A domain propagator for a constraint c is a function dom(c) mapping do- 
mains to domains which satisfies the following identity: 



Lemma 2.1. A domain propagator dom(c) for a constraint c is correct, checking, mono- 
tonic, and idempotent. 



Example 2.3. Consider the constraint v C w and domain D of Example 2.2. The domain 
propagation D' = dom(c)(D) returns domain D' where D'(v) = {{1}, {1,3}} and D'(w) = 
{{1, 2}, {1, 3}}. The missing values are those that did not take part in any solution. 

Domain consistency may be difficult to achieve for set constraints, so instead we often 
need a weaker notion of consistency. The notion of set bounds consistency is commonly 



{9 I 9 G D] n solns{c) = {9 \ 9 G /(£>)} n solns(c) 




Proof. Straightforward from the definitions. 



□ 
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used. A domain D is bounds consistent for a constraint c if for every variable v G vars(c) 
the upper bound of D{v) is the union of the values of v in all solutions of c in D, and the 
lower bound of D(v) is the intersection of the values of v in all solutions of c in D. As 
above, a set of propagators F is said to maintain set bounds consistency for a constraint c 
if solv(F,D) is bounds consistent for all domains D. 

Definition 2. A set bounds propagator for a constraint c is a function sb(c) mapping 
domains to domains satisfying the following identity: 



Lemma 2.2. A set bounds propagator sb(c) for a constraint c is correct, checking, and 
idempotent. A propagator sb(c) is also monotonic on range domains. 

Proof. All of these follow trivially from the properties of dom(c) and from the extension 



Clearly for any constraint c and any v £ vars(c) we have dom(c)(D)(v) C sb(c)(D)(v) 
for all domain propagators dom(c) and all set bounds propagators sb(c). 

2.4 Boolean Formulae and Binary Decision Diagrams 

We use Boolean formulas extensively to model sets and set relations. In particular, we make 
use of the following Boolean operations: A (conjunction), V (disjunction), -i (negation), — > 
(implication), <-> (bi-directional implication), © (exclusive OR), and 3 (existential quantifi- 
cation). We use the shorthand 3yF for 3 X1 ■ ■ ■ 3 Xn F where V = {x\, . . . ,x n }, and we use 
3yF to mean 3 V /F where V = vars(F) \ V. 

Binary Decision Trees (BDTs) are a well-known method of modeling Boolean functions 
on Boolean variables. A Binary Decision Tree is a complete binary tree, in which each 
internal node n(v, t, f) is labelled with a Boolean variable v and each leaf node is labelled 
with a truth value (false) or 1 (true). Each internal node corresponds to an if-then-else 
test of the labelled variable, and has two outgoing arcs — the "false" arc (to BDT /) and 
the "true" arc (to BDT t). In order to evaluate the function represented by the binary tree, 
the tree is traversed from the root to a leaf. On reaching any internal node, the value of 
the Boolean variable whose label appears on the node is examined, and the corresponding 
arc is followed. The traversal stops upon reaching a leaf node, whereupon the value of the 
function is taken to be the label of that node. 

A Binary Decision Diagram (BDD) is a variant of a Binary Decision Tree, formed by 
relaxing the tree structure requirement, instead representing functions as directed acyclic 
graphs. In a Binary Decision Diagram, a node is permitted to have multiple parents, as 
opposed to a tree structure which requires any node to have at most one parent. This 
effectively permits common portions of the tree to be shared between multiple branches, 
allowing a compact representation of many functions. 

The primary disadvantage of both Binary Decision Trees and Binary Decision Diagrams 
is that the representation of a given function is not canonical (i.e. one function may have 
many representations). Two additional canonicity properties allow many operations on 




and idempotence properties of conv and ran. 



□ 
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BDDs to be performed efficiently (Bryant, 1986, 1992). A BDD is said to be reduced if it 
contains no identical nodes (that is, there are no nodes with the same label and identical 
"then" and "else" arcs), and no redundant tests (there are no nodes with both "then" and 
"else" arcs to the same node). A BDD is said to be ordered if there is a total order -< 
of the variables, such that if there is an arc from a node labelled v\ to a node labelled V2 
then v\ -< V2- A Reduced Ordered BDD (ROBDD) is canonical function representation up to 
reordering, which permits an efficient implementation of many Boolean function operations. 
For more details the reader is referred to the work of Bryant (1992), or the introduction of 
Andersen (1998). 

We define the size \R\ to be the number of non-leaf nodes of the ROBDD R, as well as 
defining VAR(R) to be the set of ROBDD variables that appear as labels on the internal 
nodes of R. We shall be interested in stick ROBDDs, in which every internal node n(v, t, f) 
has exactly one of t or / as the constant node. 

Example 2.4. Figure 1(a) gives an example of a stick ROBDD LU representing the formula 
«3 A -1V4 A -1V5 A vq A v-j. \LU\ = 5 and VAR(LU) = {v^, v&, v$, vq, vj}. Figure 1(b) gives 
an example of a more complex ROBDD representing the formula —>(vi <-* vg) A —>(v2 <-> vg). 
\R\ = 9 and VAR(R) = {vi,V2,v$,vg}. One can verify that the valuation {v\ 1— > 1,^2 1— > 
0, vg 1 — >• 1 , V9 1— > 0} makes the formula true by following the path right, left, right, left from 
the root. 

2.5 ROBDD Operations 

There are efficient algorithms for many Boolean operations applied to ROBDDs. The 
complexity of the basic operations for constructing new ROBDDs is 0(|i?i||i?2|) for R± AR2, 
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node(v, t, f) if (t = /) return t else return n(v, t, f) 

and(Ri,R 2 ) 

if = or R 2 = 0) return 
if (R\ = 1) return i?2 
if (i?2 = 1) return i?i 

n(v 2 ,t 2 ,f 2 ) ■= Ri 

if (ui ^ U2) return node(vi,and(ti, i?2),and(/i, .R2)) 
else if (v 2 ^ «i) return node(v 2 ,and(t 2 , R\),and(f 2 , Ri)) 
else return node(-ui,and(ti, i2),and(/i, f 2 )) 

exists(f , R) 

if (R = 0) return 
if (R = 1) return 1 

n(v r ,t, f) := R 

if {v r ^ f) return node(f r ,exists(w, t),exists(v, /)) 
else if (v <v r ) return R 
else return or(t, /) 

Figure 2: Example ROBDD operations 

R\ V R 2 and i?i <-> i?2, 0(|i?|) for and 0(|i?| 2 ) for 3f i?. Note however that we can test 
whether two ROBDDs are identical, whether i?i <-> i?2 is equivalent to true (1), in 0(1). 

We give code for conjunction R\ A R 2 = ax\d{R\,R 2 ) and existential quantification (of 
one variable) 3v R = exists(u, R). The code for disjunction R\ V R 2 = oi(Ri,R 2 ) is dual 
to and and is very similar in structure. The code make use of the auxiliary function node 
which builds a new ROBDD node. The node function returns tift = f, and is in practice is 
memoed so any call to node with the same arguments as a previous call returns a reference 
to the previously created ROBDD node. 

Modern BDD packages provide many more operations, including specialized implemen- 
tations of some operations for improved speed. Some important operations for our purposes 
are existential quantification of multiple variables V for a formula R3V R, and the combi- 
nation of conjunction and existential quantification 3V R\ A R 2 . 

Although in theory the number of nodes in an ROBDDs can be exponential in the 
number of variables in the represented Boolean function, in practice ROBDDs are often 
very compact and computationally efficient. This is due to the fact that ROBDDs exploit 
to a very high-degree the symmetry of models of a Boolean formula. 

3. Modeling Set CSPs Using ROBDDs 

In this section we discuss how to solve set constraint satisfaction problems using ROBDDs. 
There are three parts to this — the modeling of set domains as ROBDDs, the modeling of 
set constraints as ROBDDs, and how to use both to produce a set solver. 
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3.1 Modeling Set Domains using ROBDDs 

Suppose the universe U = {1, . . . , N}. We assume that all set values are subsets of U. 
Let x be a set variable, and let D be a domain over V. Given that the size of the uni- 
verse is bounded, we can associate a Boolean variable x% with each potential element 
i G {1, . . . ,N} of x. Hence a set variable x is represented by a vector of Boolean vari- 
ables V(x) = (xi, . . . , Xjv). 

Take any set A G D(x). Then we can represent A as a valuation on the variables 
(xi,...,xjv) defined as 9 a = {x\ i— > (1 € A),...,x n t— ► (n G A)}. We can represent 
this valuation 9 a and consequently the set A as a Boolean formula B(A) which has this 
valuation as its unique solution: 



B{A) = f\ yi where y { 



ieu 



if i G A 
otherwise 



Given that we can represent any element of D(x) as a valuation, we can represent 
D(x) itself as a Boolean formula B(D(x)) whose solutions 9a correspond to the elements 
A G D(x). That is, B(D{x)) is the disjunction of B(A) over all possible sets A G D(x): 

B(D(x)) = \J B(A) where B(A) is defined above 

AeD(x) 

Any solution 6{x) G D(x) then corresponds to a satisfying assignment of the Boolean 
formula B(D{x)). From now on we will overload the notion of a domain D to equivalently 
return a set of possible sets D(x) for a variable x, or its equivalent Boolean representation 
B(D(x)). 

Example 3.1. Let U = {1, 2, 3}. Suppose v is a set variable with D(y) = {{1}, {1, 3}, {2, 3}}. 
We associate Boolean variables {^1,^2,^3} with v. We can then represent D(v) as the 
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Boolean formula (v\ A -1V2 A -1U3) V (vi A -1U2 A V3) V A v% A U3). The three solutions to 
this formula correspond to the elements of D{v). The corresponding ROBDD B(D(v)) is 
shown in Figure 3(a). 

This representation is useful since such Boolean formulae can be directly represented as 
reasonably compact ROBDDs. Given a Boolean formula representing a domain, we can 
clearly construct the corresponding ROBDD in a bottom-up fashion, but in practice we 
only ever construct the ROBDD for a domain implicitly through constraint propagation. 
As we will show, ROBDDs permit a surprisingly compact representation of (many) subsets 
ofV(U). 

Since ROBDDs are ordered, we must specify a variable ordering for the Boolean variables 
that we use. We arbitrarily order the ROBDD variables for a single set variable x as 
follows: x\ -< X2 -< ... -< x n . Assuming there is no special relationship between the 
elements of the universe, the choice of the ordering of ROBDD variables that represent a 
particular set variable is unimportant. However, the relative ordering of ROBDD variables 
comprising different set variables has a drastic effect on the size of the ROBDDs representing 
constraints, and is discussed in Section 3.2. 

Example 3.2. An ROBDD representing the 2 N ~ 5 sets in the interval [{1, 3, 4, 5}, {1, 3, 4, 5, 6, 
is shown in Figure 3(b). 

The ROBDD representation is flexible enough to be able to represent compactly a wide 
variety of domains, even those that might at first seem difficult to represent. For example, 
the set of all subsets of {1,2,3,4,5} of cardinality 2 can be represented by the ROBDD 
shown in Figure 3(c). 

For convenience, we set all of the initial variable domains Di n u(x) = V(U), which 
corresponds to the constant "1" ROBDD. Restrictions on the initial bounds of a domain 
can instead be expressed as unary constraints. 

3.2 Modeling Primitive Set Constraints Using ROBDDs 

A major benefit of using ROBDDs to model set constraint problems is that ROBDDs can 
be used to model the constraints themselves, not just the set domains. Any set constraint 
c can be converted to a Boolean formula B(c) on the Boolean variables comprising the set 
variables vars(c), such that B(c) is satisfied if and only if the corresponding set variable 
valuations satisfy c. As usual, we represent B(c) as an ROBDD. 

Example 3.3. Let v and w be set variables over the universe U = {1,2,3}, and let c 
denote the constraint » C w. Assume that the Boolean variables associated with v and 
w are (vi, V2, U3) and (w\, W2 ,103) respectively. Then c can be represented by the Boolean 
formula (v± — > w±) A (v^ — > 102) A (^3 — > W3). This formula can in turn be represented by 
either of the ROBDDs shown in Figure 4 (depending on the variable order). 

Example 3.4 demonstrates that the order of the variables within the ROBDDs has a 
very great effect on the size of the formula representations. 

Example 3.4. Consider again the constraint v C w as described in Example 3.3. Figure 4 
shows the effect of two different variable orderings on the size of the resulting BDD. In this 
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(a) Order: v\ -< V2 -< v% ~< w\ -< W2 -< W3 (b) Order: v\ -< w 1 -< V2 -< W2 ~< V3 ~< W3 

(Arcs to '0' are not shown for clarity) 



Figure 4: Two ROBDDs for v C w 



case, the variable ordering v± ~< w± -< vi -< W2 gives a much more compact representation 
of constraint than the ordering v\ -< V2 -< w± ~< W2- As can be seen from the figure, the 
latter ordering has size exponential in n, whereas the former has size linear in n. 



If the variable set V = {v±, . . . , v m } has corresponding Boolean variables (f 1,1, ^1,2, ■ • , vi,n), 
(^2,1, ^2,2, • • • , V2,n), ■ ■ ■ , {v m ,i,v mj 2, • • • , ^m,Ar), we choose to order the Boolean variables 
vi,i -<■••-< f m ,i -< ^1,2 -<■■•-< v m ^ -<■■■< v\ t N -<■■•-< v m> N- This ordering guaran- 
tees a linear representation for all of the primitive set constraints except cardinality. The 
reason is that primitive set constraints except cardinality are defined elementwise, that is 
an element i G v never constrains whether j £ w or j ^ w for i / j. For this reason 
there are no interactions between bits v\ and Wj, i 7^ j. Placing all the bits for the same 
element adjacent in the order means that the interactions of bits Vi and Wi are captured in 
a small ROBDD, which is effectively separate from the ROBDD describing the interactions 
of the next element's bits. Table 1 contains a list of the primitive set constraints, their 
corresponding Boolean formulae and the sizes of the ROBDDs representing those formula? 
under this point-wise ordering. 

Cardinality constraints can be represented in a quadratic number of ROBDD nodes 
using a simple recursive definition. We define card ((v^, ... ,Vi n ) ,l,u) to be the Boolean 
formula which restricts the number of true bits in the vector (vi 1 , . . . , Vi n ) to between I and 
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Bool pan pynrpssion F>(c\ 


size of ROBDD 


h £ v 


U K 




k v 


—<Vl. 


0(1) 


u = d 


A -^j Ui A A 1 ^ -s- AT -HJ ~<Ui 

1 \i£d 1 1 \l<i<N,i<£d 1 


O(N) 


u — V 


Ai^ -^AiiUi <— > Vi\ 
1 \l<i<N \ 1 l > 




u CI v 


1 \l<i<N\ u ' 1 Ul > 


CHAP) 


11 = 7? 1 J 7/? 


A. ,. < — > [ 7; - \/ 7 / J ■ ) 1 
1 \l<i<N\ ai \ Ul v "V/ 1 




ii = v n t(j 


Ai ^ AiiUi <— > ffi A WiS\ 
1 \l<i<N \ 1 \ % % l ) 




u = v — w 


f\l<i<N( U i (ViA^Wi)) 


0(N) 


U = V 


f\l<i<N^( u i ^ v i) 


O(N) 


U ^ V 


Vl<i<N^( u i v i) 


O(N) 


\u\ = k 


card(V(u), k, k)) 


0{k{N -k)) 


\u\ > k 


card{V{u),k,N) 


0(k(N-k)) 


\u\ < k 


card(V(u),0,k) 


0{k(N-k)) 



Table 1: Boolean representation of set constraints and the size of the corresponding 
ROBDD. 



u inclusive. 



card((vi 1} . ..,v in ),l,u) 



1 if I < A n < u 

if n< Z V u < 

(->v il A card((v i2 ,. . . , v in ) , I, u)) V otherwise 
(v h A card((v i2 ,. . .,v in ),l -l,u-l) 



It is clear from the structure of card^v^ , . . . , Vi n ) , /, u) that the resulting ROBDD is 0(n 2 ) 
in size. A more general method of characterising cardinality constraints will be presented 
in Section 5.5. 

3.3 A Basic Set Constraint Solver 

We now show how to construct a simple set domain propagator dom{c) for a constraint c. 
If vars(c) = {v±, . . . ,v n }, then we define a function dom(c) mapping domains to domains 
as follows: 

dom(c)(D)(v t ) = {M^) B ^ A A?=i DM if ^ e vars(c) 
)D(vi) otherwise 

In other words, to perform propagation we take the conjunction of the Boolean represen- 
tations of the current domains of the variables in vars(c) with the Boolean representation 
of the constraint B(c) and project the result onto the Boolean variables V(vi) representing 
each variable V{. Since B{c) and all D{vi) are ROBDDs, this formula can be implemented 
directly using ROBDD operations. 
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Figure 5: ROBDDs used in the domain propagation of c = v C w (a) D(w), (b) B{c) A 
D(v) A D(w), (c) D'(v), and (d) D'H 



Example 3.5. Consider the domain propagation of constraint v C w; with initial domain 
= {{1}, {1,3}, {2, 3}}, and D(w) = {{2}, {1, 2}, {1, 3}} from Examples 2.2 and 2.3. 
We conjoin the ROBDD B(c) shown in Figure 4(b) with the domain ROBDD D{v ) shown 
in Figure 3(a) and the ROBDD for D(w) shown in Figure 5(a). The result is an ROBDD 
representing solutions of the formula v C w A v G D(v ) Aw £ D{w) shown in Figure 5(b). 
We project the resulting onto V(v) and V(w), individually obtaining the ROBDDs D'{v) = 
{{1},{1,3}} and D'(w) = {{1, 2}, {1, 3}} shown in Figure 5(c) and (d) respectively. 

We need to verify the correctness of the propagator: 

Lemma 3.1. Let V be the collection of all set variables, and let c be a set constraint with 
vars(c) = {vi, ...,^}CV. Then dom(c) is a domain propagator for c. 

Proof. We need to verify that dom(c) satisfies the identity of Definition 1. Suppose D is a do- 
main over V. We need to check for each v G vars(c) that {char(9(v)) \ 6 G D A 6 G solns(c)} = 
{u | u \= 3y( Vi }B(c) A AILi D(vi)}, where char(X) denotes the characteristic vector of X. 
This is clearly true, since the values 9 G solns(c) are by definition the satisfying assignments 
of B(c), and the values 9 G D are by definition the satisfying assignments of A™ =1 D(fj). 
Hence the equality holds, implying that dom(c) is a domain propagator. □ 

If F is the set of all such domain propagators, then we can define a complete propaga- 
tion algorithm solv(F,D): For simplicity we use the set of constraints C rather than the 
corresponding propagators. 

solv(C, D) 
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Q ■■= C 

while (3 c£ Q) 

D' := dom(c)(D) 
if (| vars(c)\ = 1) C := C - {c} 
V := {v G V | £>(«) / £>»} 
Q:= (QU{c'GC| uors(c') D V ^ 0}) - {c} 
D := D' 
return D 

We maintain a queue Q of constraints to be (re-)propagated, initially all C. We select a 
constraint c from the queue to propagate, calculating a new domain D'. If the constraint is 
unary we remove it from C, never to be considered again, since all information is captured 
in the domain. We then determine which variables V have changed domain, and add to 
the queue constraints d € C involving these variables, with the exception of the current 
constraint c. 

By combining this algorithm with the modeling techniques described earlier, we have 
shown how to construct a simple ROBDD-based set domain propagator. Various improve- 
ments to this basic scheme will be discussed in Section 4 and Section 6. 

The reader may wonder whether in this section we have done anything more than map set 
constraints to Boolean constraints and then apply a Boolean ROBDD solver to them. This 
is not the case. Crucially the domain propagation is on the original set variables, not on the 
Boolean variables that make them up. In fact the set bounds consistency approach (Puget, 
1992; Gervet, 1997) can be considered as simply a mapping of set constraints to Boolean 
constraints. 

4. Effective Modeling of Set Constraints Using ROBDDs 

In this section we demonstrate that the ROBDD-based modeling approach is very flexible, 
allowing us to produce highly efficient implementations of many complex constraints. 

4.1 Combining Constraints and Removing Intermediate Variables 

It is rarely possible to express all of the constraints in a real set problem directly as primitive 
set constraints on the variables of the original problem. Instead, we would usually like to 
express more complicated set constraints, which can be decomposed into multiple primitive 
set constraints, often requiring the introduction of intermediate variables. 

Example 4.1. Let c be the constraint \v l~l w\ < k, which requires that v and w have at 
most k elements in common. This constraint is used in modeling the problem of finding 
Steiner systems (see Section 7.1). Since this is not a primitive set constraint, existing 
set solvers would usually implicitly decompose it into two primitive set constraints and 
introduce an intermediate variable u. The representation of the constraint then becomes 
3u u = v n w A \u\ < k. 

In the case of Example 4.1 the decomposition does not affect the strength of the resulting 
propagator. To prove this fact, we use two results presented by Choi, Lee, and Stuckey 
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(2003). Choi, Lee, and Stuckey prove these results in the case of a finite integer domain 
solver, although the proofs in the set domain case are identical. 1 

Lemma 4.1 (Choi et al., 2003). Let c\ and C2 be set constraints. Then solv{{dom(c\ A 
C2)},D) C solv({dom(ci), dom(c2)},D) for all domains D. 

Lemma 4.2 (Choi et al., 2003). Let c\ and C2 be two set constraints sharing at most 
one variable x £ V. Then solv({dom(c\), dom(c2)},D) = solv({dom{c\ A C2)},D) for all 
domains D. 

Even if the strength of the propagator is unaffected by the decomposition, splitting a 
propagator introduces a new variable, thus slowing down the propagation process. In cases 
where two set constraints share more than one variable Lemma 4.2 does not apply, so in 
such a case there can also be a loss of propagation strength. 

The ROBDD representation of the constraints allows us to utilise such complex con- 
straints directly, thus avoiding the problems associated with splitting the constraint. We 
can directly construct an ROBDD for complex constraints by forming the conjunction of 
the corresponding primitive constraints and existentially quantifying away the intermediate 
variables. 

Example 4.2. Consider the constraint c = \v PI w\ < k as discussed in Example 4.1. 
We can build a domain propagator dom(c) directly for c by constructing the ROBDD 
3V(u) u = v n w A \u\ < k. In this case the size of the resulting ROBDD is O(kN). 

The ROBDD for \v n w\ < 2, that is 3V(u) B(u = v n w) A B(\u\ < 2), where U = 
{1,2,3,4,5} is shown in Figure 6(a). The ROBDD for |m| < 2 is shown in Figure 6(b) for 
comparison. Note how each m node in the Figure 6(b) is replaced by Vi A Wi (the formula 
for Ui) in Figure 6(a). 

4.2 Modeling Global Constraints 

As the previous section demonstrates, it is possible to join the ROBDDs representing prim- 
itive set constraints into a single ROBDD representing the conjunction of the constraints. 
We can use this to join large numbers of primitive constraints to form global constraints, 
which in many cases will improve performance due to stronger propagation. 

It is one of the strengths of the ROBDD-based modeling approach that it is trivial to 
construct global constraints simply using ROBDD operations on primitive set constraints, 
without laboriously writing code to perform propagation on the global constraint. This 
approach is very powerful, but is not feasible for all combinations of primitive constraints. 
As we shall see, some global constraints that we might desire to construct lead to ROBDDs 
that are exponentially large. 

A useful global constraint is the constraint partition(t> i, . . . , v n ), which requires that 
the sets v\,V2, ■ ■ ■ ,v n form a partition of the universe U = {1, . . . , iV}. We can easily 



1. The observation of Choi et al. that Lemma 4.2 does not apply if the shared variable is a set variable is 
only true if we are performing set bounds propagation, not set domain propagation. 
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construct this constraint from primitive set constraints as follows: 

n— 1 n 

partition^!, . . . ,v n ) = f\ /\ 3 Uij (uij = V; t n Vj A = <p) 

i=l j=i+l 

n 

A 3 W() ■ ■ ■ 3 Wn (w = 4> A (/\ Wi = wi-i U Vi) Aw n =U) 

i=\ 

The propagator dom(partition(- • • )) is stronger than a domain propagator on the 
decomposition. For example, consider the constraint c = partition(x, y, z) as depicted 
in Figure 7, and the domain D where D(x) = D{y) = {{1},{2}} and D(z) = {{2}, {3}}. 
Domain propagation using the decomposition of the constraint will not alter D, but domain 
propagation using the global constraint will give dom(c)(D)(z) = {{3}}. Hence stronger 
propagation can be gained from global propagation using this constraint. 

Unfortunately not all global constraints can be modelled efficiently using this approach. 
In particular, there is a risk that the ROBDD representation of a global constraint could 
be extremely large, making it infeasible to construct and use for propagation. 

For example, consider the constraint atmost^Vi, . . . , v n ) , k) proposed by Sadler and 
Gervet (2001), which requires that each of the sets v\,...,v n has cardinality k and the 
intersection of any pair of the sets is at most 1 element in size. This constraint models the 
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Figure 7: ROBDDs for the constraints partition(x, y, z) and lexlt(x,y), where x, y and z 
are set variables taking values from the universe U = {1, 2, 3}. 



Steiner triple systems of Section 7.1. Bessiere, Hebrard, Hnich, and Walsh (2004) proved 
that enforcing bounds consistency on this constraint is NP-hard, so it follows that enforcing 
domain consistency on this constraint is at least NP-hard as well. Theoretically we can still 
construct an ROBDD representing this constraint from primitive constraints as follows: 

n n—1 n 

f\ \vi\ = n A /\ /\ 3 UlJ (uij =Vif\ vj A \uij\ < 1) 

i=l i=l j=i+l 

Unfortunately, the resulting ROBDD turns out to be exponential in size, making it im- 
practical for use as a global propagator. This is not surprising in light of NP-hardness of 
the problem in general — in fact it would be surprising if the resulting ROBDD was not 
exponential in size! 

4.3 Avoiding Symmetry Through Ordering Constraints 

It is important in modeling a constraint satisfaction problem to minimize symmetries in 
the model of the problem. A model that contains symmetrical solutions often has a greatly 
enlarged search space, leading to large amounts of time being spent in searching sections 
of a search tree that are identical up to symmetric rearrangement. It is therefore highly 
desirable to remove whatever symmetries exist in a problem. 
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One approach to symmetry-breaking is the introduction of additional ordering con- 
straints between the variables of the problem. A convenient ordering to use on sets is a 
lexicographic order on the characteristic bit vectors of the sets. In other words, if v and w 
are set variables, then v < w in the lexicographic ordering if and only if the list of bits V(v) 
is lexicographically smaller than V(w). We can model this lexicographic ordering constraint 
as lexlt(v,w, 1), which is defined recursively in the following manner: 



An example of such an ROBDD is depicted in Figure 7. 

We shall make use of the lexicographic ordering constraint extensively in our experiments 
in Section 7. 

5. Modeling Integers, Multisets, and Weighted Sum Constraints 

In this section we show how to model integer variables and integer constraints using ROB- 
DDs. In general such representations can be very large, which limits the usefulness of 
ROBDDs as the basis for a general purpose finite-domain constraint solver. Despite this, 
the ability to represent integers is extremely useful as a component of a set solver. We 
propose two major uses of the integer representation. 

Firstly, we can use the integer representation to model values such as the weighted sum 
of the elements of a set. Constraints on the weighted sum of elements of a set have been 
shown to be useful in practical applications (Mailharro, 1998). 

Secondly, we can model finite multisets using ROBDDs by replacing the individual 
ROBDD variables of the set representation with bundles of ROBDD variables, each bun- 
dle corresponding to a binary integer. Multiset operations can then be constructed by 
composing integer operations on the variable bundles. 

In addition, the integer representation as described here could be used as an inter- 
face between an ROBDD-based set solver and a more conventional integer finite-domain 
solver. Such an interface could easily be implemented using channeling constraints between 
ROBDD integer and finite-domain versions of the same variable. 

5.1 Representing Integer Values using ROBDDs 

In order to model integers and integer operations we must choose an appropriate represen- 
tation in terms of Boolean formulas. In general we are free to use any encoding of an integer 
as a binary sequence, such as unary, unsigned binary and twos-complement encodings, but 
for simplicity we choose to represent all of our integers in unsigned binary form. 

We can represent an arbitrary integer expression e by a list of Boolean formulas. Each 
formula in the list corresponds to a single bit of the unsigned binary value of the expression. 
We will denote such a list by (e n _i, e n _2, . . . , ei, eo), where each is a Boolean formula, in 
order from the most significant bit to the least significant bit. We interpret a formula as 
a "1" bit if the formula is logically true, and as a "0" bit otherwise; for simplicity we will 
just call the formulas the "bits of the expression". We will also use the expression e and its 
list of constituent bits interchangeably. As usual, we will represent the Boolean formulas 




if n > N 

(->v n A w n ) V {(v n <-> w n ) A lexlt(v, w, n + 1) otherwise 
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as ROBDDs. As we shall see, this notation is flexible enough to represent arbitrary integer 
expressions. 

Example 5.1. Consider the integer constant k = 25, which is 11001 in unsigned binary. 
We will represent k as the list (1,1,0,0,1). 

In order to represent an integer variable x, we associate with x a fixed set of Boolean 
variables {x^-i, . . . , xo}. The value of x is then taken to be the value of the unsigned binary 
integer (x^-i, ■ ■ ■ , xo). By varying the value of the Xj variables, the value of x can range 
from to 2 k — 1 inclusive. As always the ordering of the Boolean variables has a significant 
effect on the sizes of the ROBDD representations of formulas. If x = (xk-i, ■ ■ ■ ,xq), y = 
{yk-i, ■ ■ ■ , yo), and z = (zk-i, • • • , zq), then we choose to order the corresponding ROBDD 
variables in a interleaved most-significant-bit-first order, i.e. Xk-i -< yu-i ~< z k-i ~< %k-2 < 
■ ■ ■ -< xo -< y -< z . 

Since we permit arbitrary Boolean formulae as the bits of an expression, we can also 
model arbitrary integer expressions. For example, suppose x and y are integer variables with 
bits (x2, x±, xo) and (j/2, yi,yo) respectively. Then, for example, we can represent the expres- 
sion xAy (where the A denotes the bitwise AND operator) as the list (x2 A 7/2, X\ A y±, xo A yo) 
We are not limited to logical operations, as we shall see in the next section. 

5.2 Representing Integer Operations using ROBDDs 

We can also use ROBDDs to model arithmetic operations such as addition, by analogy 
with the design of the corresponding logic circuits. For the purpose of the set and multiset 
solvers, we only require implementations of the operations of addition, left shift, minimum, 
maximum, multiplication by a constant, and multiplication by a single variable bit. We do 
not require a general implementation of multiplication using ROBDDs. 

It is convenient to assume that all integer expressions have the same number of bits. We 
may assume this without loss of generality since we can freely pad the left of the shorter of 
any pair of expressions with "0" bits. 

To model addition, we simulate the operation of a full binary adder. Suppose x and y 
are integer expressions, with bit representations (x;_i, . . . , xo) and {yi-i, . ■ ■ , yo}- We can 
use ROBDDs to compute the output bits plus(x, y) of the operation x + y as follows (here 
Cj denotes a carry bit and Si denotes a sum bit): 

c_i = 

Si = Xi (B yi (B Cj_i 

Ci = (->Cj_i A Xi A yi) V (cj_i A (x, V y^) 
plus(x,y) = (q_i,sj_i,...,si,s ) 

Note we avoid overflow by extending the size of the result by one bit. 

Example 5.2. Suppose x is an integer variable, with bits (xi,xo). We can represent the 
expression x+3 by the bits ((^xo A x±) V xo), x\ © 1 © xo, xo © 1) = (xo V xi, xo «-» x\, -ixo). 

The left shift operation is trivial to implement. If x is an integer expression represented 
by (x/_i, . . . , Xq) and k is a non-negative integer, then we can represent the left shift of x 



< i < I 
0<i<l 
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Figure 8: ROBDDs representing the bits (e3, e2, ei, eo) of the expression e 
mul((xi,xo) ,3), together with the simplified Boolean expressions for each e$ 



by /c bits sM(x, fe) by the following formula: 
shl((xi-i, ...,x ),k)-- 



(xi- 



k bits 



We can implement the operation of multiplication by a constant using the plus and shl 
operators. If x is an integer expression with bits (xi—i, . . . , xq) and k is a non-negative 
integer, then x x k corresponds to mul(x, k) in the following formula: 



mul(x, k) 



r<> 

mul(shl(x, 1), |) 
^plus(x, mul(shl(x, 1), |_f J)) if & is °dd and fe > 



if fe = 

if fe is even and fe > 



(1) 



Example 5.3. Let x = (xi,xq), and consider the expression e = mul(x,3). By applying 
Equation (1), we obtain e = (xi A (xi A xo), x\ © (xo A x±), xq ffi x\, xq) which can be sim- 
plified to (xo A x±, -1X0 A xi, xo ffi xi, xo). The corresponding ROBDD representations are 
shown in Figure 8. 



5.3 Integer Constraints using ROBDDs 

We can also express constraints on integer expressions using ROBDDs. In particular, we 
show how to implement equality and inequality constraints. As usual, we assume any two 
expressions x and y have equal lengths of I bits; if not, we pad the shorter expression with 
"0" bits on the right. 

Equality of two integer expressions x and y is easy to represent as an ROBDD — the 
corresponding bits of x and y must be equal. Hence we can represent the equality constraint 
x = y as the ROBDD B(x = y): 

i-i 

B(x = y) = f\xi<r+yi 

i=0 
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Note that this is identical to the implementation of equality for two set expressions with 
the addition of zero-padding. 

It turns out that we already have an implementation for inequality constraints, albeit 
in a different guise. Inequalities of the binary integers correspond to inequalities on the 
lexicographic ordering of the bit representations, so we can implement, for example, the 
strict less-than constraint x < y for two integer variables x and y using the lexlt operation 
from Section 4.3. 

B{x < y) = lexlt(x,y) 

We can then use this to construct the reverse inequality by swapping the order of the 
operands, and the non-strict inequality by negating the formula (reversed as necessary). 

The implementation of inequalities also leads us to an implementation of the minimum 
and maximum expressions. Consider the problem of finding the smaller of two integer 
expressions x and y. If x and y have bit vectors (xi-i, . . . , xq) and (yi-i, ■ ■ ■ , yo) respectively, 
we recursively define min(x, y) as follows: 

U = Ri = 

Li = Lj+i V (->L i+ i A -i-Rj+i A -iXj A yi) 1 < i < I 

Ri = R i+1 V (->L i+ i A -^R i+ i A Xi A ->?/;) 1 < i < I 

nii = A Xi) V (Ri+i A yi) V (->Lj+i A ^R i+ i A Xi A yi) <i <l 

min(x,y) = (mj_i, . . . ,m ) 

In the equation above, the Li and Ri values are flag bits which state whether the higher 
order bits have already allowed us to conclude that one of the two values is the minimum. 
The maximum operation is defined similarly. 



5.4 Modeling Multisets and Multiset Constraints 

Various authors have suggested that multisets are a valuable addition to the modeling 
abilities of a set constraint solver (Kiziltan &; Walsh, 2002). In this section we briefly show 
how multisets and multiset constraints can be modelled using ROBDDs by making use of 
the integer building blocks described above. 

A multiset m is a unordered list of elements {mo, . . . , m n ^ drawn from the universe U, 
in which (unlike a set) repetition of elements is permitted. Most set operations have parallel 
operations on multisets, although the multiset operations are not strict generalisations of 
the set operations. Let occ(i, m) denote the number of occurrences of an element i in a 
multiset m. Suppose m and n are multisets, and k is an integer constant. We define the 
following multiset relations and operations by their actions on the number of occurrences 
of each element i in the universe: 

• Equality: m = n iff occ(i, m) = occ(i, n) for all i &U. 

• Subset: m C n iff occ(i, m) < occ(i, n) for all i €U. 

• Union: occ(i, m U n) = occ(i, m) + occ(i, n) for all i GW. 

• Intersection: occ(i, mDn) = mm{occ(i, m), occ(i, n)} for all i GU. 
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• Difference: occ(i, m\n) = max{0, occ(i, m) — occ(i, n)} for all i & 

• Cardinality: |m| = J2i<^u °cc(i,m) for all i € U. 

To represent a set variable x we associated a vector of Boolean variables {x±, . . . ,x n ) 
with the bits of the characteristic vector of a valuation of x. In the case of a multiset, the 
characteristic vector of a multiset m is a vector of integers, and so we need to associate an 
integer value mi with each potential element i of the multiset m. We can model each such 
integer value using the approach described above. 

If m is a multiset variable, then we associate a bundle of ROBDD variables with every 
i £U, the contents of which comprise the bits of the corresponding integer expression wij. In 
order that we can represent multisets in a finite number of bits, we assume that the number 
of occurrences of any element in a multiset variable is bounded above by a reasonably small 
M, allowing us to use only k = |"log 2 M] Boolean variables per bundle. We then write m 
as a list of bundles (mi, . . . , m n ), where each mi is in turn a list of bits (m^k-i, • • • , rriifi}- 

Given a representation of multiset variables, we now turn our attention to the implemen- 
tation of multiset expressions and constraints. Multiset expressions can be implemented in 
the obvious way as sequences of integer expressions. For example, suppose x and y are mul- 
tiset variables with associated bundles (x±, . . . ,xn) and {y±, . . . , yw) respectively. Then the 
bundles corresponding to the expression xUy are (plus(xi, yi), plus(x2, 2/2), • • ■ , plus(xN, Un)}- 
Similarly, the bundles corresponding to the expression xf]y are (min(xi, y±), . . . , min(a;jv, Un)), 
and so on for other expressions. We show how to implement cardinality and weighted sum 
constraints in Section 5.5. 

Multiset constraints are also trivial to implement — for example, two multisets x and y 
are equal if and only if all of their constituent bundles are equal, and so multiset equality 
can be modelled by a conjunction of integer equalities. Relations such as subset correspond 
to a conjunction of integer inequalities on the constituent bundles, the implementation of 
which was described in Section 5.3. 

Up until this point we have left the ordering of the ROBDD variables that make up a 
multiset variable unspecified. Unfortunately, unlike the set case, there is no single optimal 
variable ordering that is guaranteed to produce compact descriptions of all the primitive 
multiset constraints. For example, a subset constraint can be compactly represented by 
a bundle-major bit ordering but not a bit-major ordering (see Figure 9), since a subset 
constraint consists of a series of integer inequalities between the corresponding bundles of 
two multiset variables, and so a bundle-major ordering gives an interleaving of variables as 
above. However, the opposite is true for a cardinality constraint, which consists of a sum 
of the values of the bundles within a variable, the variables of which are interleaved under 
a bit-major ordering but not a bundle-major ordering. These two orderings are mutually 
exclusive, and hence we can conclude that in general there need not be an optimal variable 
ordering for modeling multiset constraints. 

5.5 Weighted Sum and Cardinality Constraints 

In many practical applications we are interested in placing constraints on a weighted sum of 
the elements of a set variable. For example, in the Balanced Academic Curriculum problem 
(problem prot>030 of CSPLib), every course has an associated weight corresponding to its 



130 



Solving Set Constraint Satisfaction Problems using ROBDDs 




Bit major Bundle major 



Figure 9: Bit and bundle major orderings for two multiset variables 



academic load, and there is a limit to the total academic load that can be undertaken during 
any given time period. In the set model of the problem proposed by Hnich, Kiziltan, and 
Walsh (2002), the limits on the academic load in a period are made through the use of 
a weighted sum constraint. In addition the cardinality constraint which we have already 
described is a special case of the weighted sum constraint where all the weights are set to 
"1" . It therefore seems essential to implement this constraint in the ROBDD framework. 

Suppose x is a set or multiset expression with bit bundles (xi, . . . ,x n ). (If x is a set 
expression, then each bit bundle has size 1). We can use the integer operations described 
earlier to produce an integer expression wsum(x, w) corresponding to the weighted sum 
Y^i=i XiWi, where it; is a vector of integers (wi, . . . , w n ): 

<Po = 

4>i = plus((/>i-i,mul(xi,Wi)) 
wsum(a;, w) = 4> n 

Expressions involving the cardinality of set and multiset variables can then be expressed 
as a special case of the weighted sum expression where Wi = 1 for all 1 < i < n. We have 
already seen one method of constructing ROBDDs for the constraints |x| = k, \x\ < k and 
\x\ > k in the set variable case; this method is more general, since it permits us to directly 
model constraints such as \x\ + 5 < |y| + \ z\. This is of great practical value — for example it 
is needed to implement the Hamming Code experiment of Section 7.3. In the cases already 
discussed in Section 3.2 the ROBDDs produced by the two methods are identical since 
ROBDDs are a canonical representation. 

6. Efficient Constraint Propagation Using ROBDDs 

In this section we discuss improvements to and variants on the basic domain propagation 
scheme presented in Section 3.3. 

The implementation of the domain propagator dom(c) has a substantial effect on the 
performance of the solver. The definition of dom(c) was given in its purest mathematical 
form for simplicity. Section 6.1 discusses several implementation details that can lead to 
greatly improved efficiency when performing domain propagation. 

In general, the inferences that can be obtained using the domain propagator may be 
too costly to be practical, and so in some circumstances it may be desirable to enforce a 
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weaker form of consistency. Weaker consistency may be substantially cheaper to enforce, 
permitting more time to be spent in searching the solution space. 

As always there is a compromise between propagation time and search time, and for 
certain problems it may be more productive to spend more time searching than performing 
more accurate propagation, while for other problems the converse may be true. The ROBDD 
based representation allows this to be taken to extremes — in theory it is possible to form 
a single ROBDD representing the solutions to a constraint satisfaction problem by forming 
the ROBDD conjunction of the constraints. A solution could then be trivially read from 
the ROBDD as a satisfying assignment. Usually such an ROBDD would be prohibitively 
expensive to construct in both time and space, forcing us to maintain less strict consistency 
levels. 

Accordingly, we show how to implement weaker levels of consistency using the ROBDD 
representation by combining the domain propagator with an approximation operation. The 
approximation operation simplifies the ROBDD representing a domain to its bounds closure 
under a suitable definition of bounds. The ROBDD representing the bounds of a domain is 
almost always smaller than the domain itself, leading to better performance for all future 
operations on that domain. As we shall see, this can lead to substantial performance 
improvements for the overall solver. 



6.1 Domain Propagation 

Let c be a constraint, with vars(c) = {v\, . . . , v n }. Section 3.3 gave the following definition 
of a domain propagator: 

n 

dom(c)(D)(v t ) = 3 v(Vl) (B(c) A /\ D{v 3 )) (2) 

i=i 

Since B(c) and D(vj) are ROBDDs, we can directly implement Equation (2) using 
ROBDD operations. In practice it is more efficient to perform the existential quantification 
as early as possible to limit the size of the intermediate ROBDDs. Here we make use of 
the efficient combined conjunction and existential quantification operation, which we'll call 
and- abstraction, provided by most ROBDD packages. 

This leads to the following implementation: 

$ = B{c) 

jJ _ )^V{vj){ D {vj) A (fP' 1 ) 1 < i,j < n, j (and-abstraction) 

[<i>i i = 3 

dom(c)(D)(v l )=D(v i )A<ft? 

The worst case complexity is still 0(\B(c) \ x YYj=i \D{vj)\) for each Vj. Clearly some of the 

computation can be shared between propagation of c for different variables since (ft? = <$?■, 
when j < i and j < i! . Even with this improvement the algorithm of Equation (3) uses 
0(n 2 ) and-abstraction operations (which experimentally have been shown to occupy the 
majority of the execution time of the set solver). 

The domain propagator implementation of Equation (3) can be significantly improved 
by observing that in the case of an n- variable constraint (n > 3) many similar sub-formulas 
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dom_divide_conquer(Z), <p, V): 
if (\V\ = 0) return D 
else if (V = {vi}) 

D(vi) := D(vi) A (j) 

return D 
else 

{vi, . . . ,v k } := V 
h := L|J 

R ■= ^V^Divi) A 3 V ( V2 )D(v 2 ) A • • • 3 V ( Vh )D(v h ) A 

L := 3 V(v h+1 )D(v h+ i) A 3y (%+2) D(^ +2 ) A • • • 3 y(l)fe)J D(?; fc ) A <f> 

D\ := dom_divide_conquer(Z), L, {v\, . . . , Vh}) 

D2 := dom_divide_conquer(L'i, R, {vh+i, ■ ■ ■ ,t>fc}) 

return D2 

Figure 10: A divide and conquer algorithm for domain propagation 



are computed. Due to the need to perform the existential quantification operations as early 
as possible, we do not have complete freedom to rearrange the order of evaluation as we see 
fit. However, a simple divide-and-conquer strategy for calculating dom(c)(D)(vi) allows us 
to perform domain propagation using just 0(n log n) and-abstraction operations. We define 
dom(c)(D) = dom_divide_conquer(D, B(c), vars(c)), where dom_divide_conquer is defined in 
Figure 10. 

6.2 Set Bounds Propagation 

As domain propagation may be prohibitively expensive to enforce for some problems, it is 
useful to investigate less strict notions of consistency. In such cases, we can speed prop- 
agation by simplifying the domains through approximation. Since set bounds under the 
subset partial ordering relation are one of the most commonly used approximations to a set 
domain, it seems natural to implement a set bounds propagator in the ROBDD framework. 
Only relatively minor changes are needed to the domain propagator to turn it into a set 
bounds propagator. 

Given the ROBDD representation of a set domain, we can easily identify the corre- 
sponding set bounds. In an ROBDD-based set domain representation, the set bounds on a 
domain correspond to the fixed variables of the ROBDD representing the domain. We say 
an ROBDD variable v is fixed if either for all nodes n(v, t, e) t is the constant node, or 
for all nodes n(v, t, e) e is the constant node, and such a node appears in every path from 
the root of the diagram to the "1" node. 

Such nodes can be identified in a single pass over the domain ROBDD, in time propor- 
tional to its size. If 4> is an ROBDD, we will write [0] to denote the ROBDD representing 
the conjunction of the fixed variables of 4>. If 4> represents a set of sets S, then [0] represents 
conv(S). An ROBDD <f> where <p = is a stick ROBDD by definition. 



133 



Hawkins, Lagoon, & Stuckey 



Example 6.1. Let <j> be the ROBDD depicted in Figure 1(c). Then [</>] is the ROBDD of 
Figure 1(a). 

Using this operation we can convert our domain propagator into a set bounds propagator 
by discarding all of the non-fixed variables from the domain ROBDDs after each propagation 
step. Suppose that D(v) is a stick ROBDD for each v £ V. If c is a constraint, with 
vars(c) = {vi, . . . ,v n }, we define a set bounds propagator sb(c) thus: 



Despite only relatively minor differences between the set bounds propagator and the 
domain propagator, the set bounds propagator is usually significantly faster than a do- 
main propagator for two reasons. Firstly, as the domains D(v ) are all sticks, all of the 
ROBDD operations are cheap, compared to operations on the possibly very large ROBDDs 
representing arbitrary domains. The entire propagator can be implemented with 0(\B(c)\) 
complexity, since all of the other ROBDDs are sticks. Secondly, we can use the updated set 
bounds to simplify the propagator ROBDD B(c). Since domains are monotonic decreasing 
in size, fixed variables will remain fixed up to backtracking, and so we can project them 
out of -B(c), thus reducing the size of the propagator ROBDD in future propagation steps. 
This leads us to the following implementation of the propagator: 



After this propagation step we can replace the representation of the constraint B(c) by 4> n 
since the fixed variables will no longer have any new impact. 

Example 6.2. Consider bounds propagation for the constraint c = v C w where N = 3. 
The ROBDD representation B(c) is given in Figure 4. Assume the domains of v and w are 
respectively [{1}, {1, 2, 3}], represented by the formula v±, and [0, {1, 2}], represented by the 
formula ^w 3 . The ROBDD <p n = 3vi3w 3 B(c) A v\ A ~~ 1103 is shown in Figure 11(a). We 
have that [</>„] = w\ A -1U3. We can project out the fixed variables v \ , w\ , V3 , W3 from B(c) 
to get a new simplified form of the constraint t>2 — ► W2 shown in Figure 11(b). 

This set bounds solver retains all of the modeling advantages of the domain solver, 
including the ability to easily conjoin and existentially quantify constraints, to remove in- 
termediate variables and to form global constraints. In some cases this permits a substantial 
performance improvement over more traditional set bounds solvers. 

Experimentally it appears that a direct implementation of Equation (4), written to 
use a divide and conquer approach to calculate the and-abstractions, is faster than an 
implementation of Equation (5), even if a divide and conquer approach is used in calculating 
the existential quantification. The former approach calculates fewer intermediate results, 
which leads to a faster propagator overall. Experimental results for the bounds propagator 
are given in Section 7. 



n 




(4) 



00 = B{c) 

4>j = ^VAR{D{vj))D{vj) A (f>j-i l<j<n 
sb(c) (D) ( Vi ) = D( Vi )A 3 v{vi) \<t> n \ l<i< n 



(5) 
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Figure 11: Set bounds propagation on the constraint v C w showing (a) resulting ROBDD 
after conjunction with domains, and (b) simplified constraint ROBDD after 
removing fixed variables. 



6.3 Split Domain Propagation 

We can combine the set bounds propagator with the domain propagator to produce a more 
space efficient split domain propagator. By separating the domain representation into fixed 
and unfixed parts, we can reduce the total size of the representation, also hopefully speeding 
propagation. 

One of the unfortunate characteristics of ROBDDs is that the size of the ROBDD 
representing a domain is highly dependent on the variable ordering. Consider an ROBDD 
representing a set domain which contains several fixed variables. If these variables do 
not appear at the beginning of the variable ordering, then the ROBDD representing the 
domain will in effect contain several copies of the sticks representing the fixed variables. 
For example, Figure 1(c) contains several copies of the stick in Figure 1(a). Since many of 
our ROBDD operations take time proportional to the product of the number of ROBDD 
nodes of their arguments, this overly large representation has a performance cost. We can 
solve this problem in two ways — either by reordering the ROBDD variables or by splitting 
up the domain representation. 

Variable reordering is capable of eliminating redundancy in the representation of any 
individual domain, but in general cannot eliminate redundancy across a set of domains. By 
reordering the ROBDD variables, we can reduce the size of a domain by placing the fixed 
variables at the beginning of the variable order, thus removing the unnecessary duplication 
for that domain. Unfortunately, the variable order is a global property of all ROBDDs in 
existence, whereas the fixed variables of a domain are a local property specific to a particular 
domain, so there may not be a variable ordering that is optimal for all of the domains in a 
problem. 

In the context of applying ROBDDs to the groundness analysis of logic programs, Bag- 
nara (1996) demonstrated that the performance of an ROBDD-based program analyzer 
could be improved by splitting up ROBDDs into their fixed and non-fixed parts. We can 
apply the same technique here. 
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We split the ROBDD representing a domain D(v) into a pair of ROBDDs (LU, R). LU is 
a stick ROBDD representing the lower and upper set bounds on D(v), and R is a remainder 
ROBDD representing the information on the unfixed part of the domain. Logically D = 
LU A R. We will write LU(D(v)) and R(D(v)) to denote the LU and R parts of D(v) 
respectively. 

The following results provide an upper bound of the size of the split domain represen- 
tation: 

Lemma 6.1. Let G be an ROBDD, and let v be a fixed variable of G. Then \3 V G\ < \G\. 

Proof. Since v is a fixed variable, either for every node n(v,t,f) in G t is the constant 
node, or for every such node / is the constant node. Since a node n(v, t, f) corresponds to 
the proposition (v At) V A/), it is clear that 3 v n(v,t, f) corresponds to simply i V/, and 
moreover since v is a fixed node one of t or / is zero. Hence the existential quantification 
of a fixed variable v simply removes all nodes labelled v from D. Since there is at least one 
such node, the result follows. □ 

Lemma 6.2. Let G be an ROBDD, LU = [G], and R = 3 VAR{W) G. Then G <-> LU A R 
and \LU\ + \R\ < \G\. 

Proof. The result G LU A R is straightforward, so we only prove the result on sizes. 

Suppose that VAR(LU) = {v±, . . . ,v n } is the set of fixed variables of G. Then, since 
|3j;iCr| < \G\, 1 + \3 V1 G\ < \G\. By repeating this operation for each Vi, we obtain n + 
\~^var(lu)G\ ^ 1^1- But as LU is a stick, trivially \LU\ = n, and as R = ^var{LU)^ ^ 
definition this is the required inequality. □ 

Note that \D\ can be 0(\LU\ x \R\). For example, considering the ROBDDs in Figure 1 
where LU is shown in (a), R in (b), and D = LU A R in (c), we have that \LU\ = 5 and 
\R\ = 9 but \D\ = 9 + 4 x 5 = 29. 

We now show how to construct a propagator on split domains. Firstly, we eliminate 
any fixed variables (as in the bounds propagator) and then apply domain propagation on 
the remainder domains R. The propagator produces a new pair (LU, R) consisting of new 
fixed variables and a new remainder. This process is shown below: 

<Po = B{c) 

fa = 3 VAR(LU(D(v,)))( LU ( D ( v j)) A 

n 

6i = 3 v{vi) (<t> n A f\R(D( Vj ))) (6) 

3=1 

Pi = LU(D( Vi )) A [Sij 
dom{c)(D)(vi) = {/3i,3 VAR{/3i) Si) 

For efficiency the Si components can be calculated using the divide-and-conquer ap- 
proach described for the domain propagator. 

The split domain representation has three main advantages. Proposition 6.2 tells us that 
the split domain representation is no larger than the original domain representation. How- 
ever, often the split representation is substantially smaller, which can lead to improvements 



136 



Solving Set Constraint Satisfaction Problems using ROBDDs 



{1,2,3,4) 

Upper bound 
(1,2,3} (1,2,4) (1,3,4) (2,3,4) 




(1,2) (1,3) (1,4) (2,3) (2,4) (3,4) 



Lower bound 



(1} (2) {3} (4) 




(} 



Figure 12: The set interval [0, {1, 2, 3, 4}] with upper and lower cardinality bounds u = 3 
and I = 2 respectively 



in propagation performance. The split solver can also use the propagator simplification 
technique from the bounds solver by abstracting the fixed variables out of the propaga- 
tor ROBDDs. Finally using the split solver we can mix the usage of domain and bounds 
propagators in the same problem. 

Experimental results for the split domain propagator are given in Section 7. 

6.4 Cardinality Bounds Propagation 

Given that we are able to model a set bounds propagator using ROBDDs, it is also appropri- 
ate to consider how we might model other levels of consistency for set constraint problems. 
One level of consistency that is commonly used (Azevedo, 2002; Miiller, 2001) is that of 
combined set bounds and cardinality consistency, in which upper and lower bounds on the 
cardinality of each domain are maintained in addition to the bounds under the subset par- 
tial ordering. This hybrid approach allows a more accurate representation of some domains, 
particularly those with constrained cardinality, which are common in set problems. 

Example 6.3. Figure 12 depicts the set interval [0, {1, 2, 3, 4}], together with lower and 
upper cardinality bounds 2 and 3 respectively. In general an interval consists of a large 
number of sets, making it a crude approximation to a set domain. Cardinality bounds 
permit a more fine grained representation by effectively allowing us to select a subset of the 
rows of the lattice diagram. 

Just as we were able to implement a set bounds propagator using the domain propa- 
gator together with a function which extracts the set bounds of a domain, so too can we 
create a combined set bounds and set cardinality propagator. Here we will extend the split 
domain solver by simplifying the "remainder" component of a split domain to an ROBDD 
representing its cardinality bounds. 
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bdd_count_cardinality(D, V): 
if (D = 0) return (oo, — oo) 
else if (D = 1) 

if (\V\ = 0) return (0,0) 

else 

{v!,v 2 , ■ ,v n ) := V 
return (0, n) 

else 

if (|V| = 0) error 
n(v, t, e) := D 
(v 1 ,v 2 ,...,v n ) := V 
if {v\ >~ v) error 
else if (v = v\) 

(k,ut) := bdd_count_cardinality(t, (v 2 , ■ ■ ■ , v n )) 

(l e ,u e ) '■= bdd_count_cardinality(e, (v 2 , ■ ■ ■ , v n )) 

return (min(Z t + 1, l e ), max(u t + 1, u e )) 
else 

(l,u) := bdd_count_cardinality(D, (v 2 , . . . ,v n )) 
return (I, u + 1) 

bdd_card_bounds(D, V): 

(l,u) := bdd_count_cardinality(D, V) 
if (I = oo or u = — oo) return 
return card(V,l,u) 

Figure 13: An algorithm to determine the cardinality bounds on the domain of a set variable 
represented by an ROBDD D, where V is the vector of bits in the Boolean 
representation of the set variable 



As before, we need a method for extracting an ROBDD representing the cardinality 
bounds of an arbitrary domain ROBDD. We perform this operation in two stages. Firstly 
we define a function bdd_count_cardinality which takes an ROBDD representing a set domain 
and returns upper and lower bounds on its cardinality. We can represent these bounds 
in ROBDD form by constructing a new cardinality constraint ROBDD as described in 
Section 3.2. 

An implementation of bdd_count_cardinality is shown in Figure 13. This function can 
be implemented to run in 0(|-D[|V|) time if dynamic programming/ caching is used to save 
the results of the intermediate recursive calls. In practice since D and V are highly in- 
terrelated, it is 0{\D\). In our implementation we utilise the global cache mechanism of 
the ROBDD library, which also permits caching of partial results between multiple calls to 
bdd_count_cardinality. 
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Figure 14: Cardinality propagation example showing (a) resulting ROBDD projected onto 
x, and (b) calculation of cardinality bounds 



We can now use this function to construct a function bdcLcarcLbounds which takes an 
ROBDD D and a list of Boolean variables V and returns a new ROBDD representing the 
bounds on the cardinality of the solutions of D. The implementation of bdcLcarcLbounds is 
shown in Figure 13. 

The algorithm for a split set bounds and set cardinality propagator for a constraint c is 
given by sbc(c) in the following equation: 

00 = B(c) 

fa = 3 VAR(LU(D(v,)))( LU ( D ( v j)) A <A?-l) 
n 

S i = 3 v{vi) (4>nA/\R(D(v j ))) (7) 

3=1 

A = LU{D( Vi )) A [Sij 
sbc{c){D)(v t ) = (J3 i ,bd<l-<ardJBOut\ds(3 VARm 6 i ,V(vi) \ VAR^))) 

Note that we only keep the cardinality bounds on the remaining non-fixed Boolean 
variables for a set variable v rather than on all the original variables V(v), since we do 
not need to consider fixed variables again, and it leads to a slightly smaller cardinality 
ROBDD. As usual Equation 7 should be implemented using a divide- and-conquer approach 
for efficiency. Experimental results for this propagator are shown in Section 7. 

Example 6.4. We illustrate set bounds and cardinality propagation on the constraint 
c = lexlt(x,y) whose ROBDD B(c) is shown in Figure 7(b). Assume the original domains 
for x and y are universal, so D(x) = D(y) = [0, {1,2,3}], represented by the ROBDD 1. 
The ROBDD for S x = 3 v{x) B(c) is shown in Figure 14(a), and f3 x = [<5J = 1. The tree 
of calculations for bdd_card_bounds(<5 x , (x±, X2, £3}) is shown in Figure 14(b). Overall the 
cardinality of x is determined to be in the range [0, 2] . 
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6.5 Lexicographic Bounds Propagation 

An alternative form of set consistency proposed by Sadler and Gervet (2004) is to maintain 
bounds under a lexicographic ordering in addition to set bounds. The lexicographic ordering 
is a total ordering on sets which embeds the subset partial ordering. Bounds under the 
lexicographic ordering alone are not sufficient to express the effects of many constraints (in 
particular the inclusion of a single element), so Sadler and Gervet constructed a hybrid solver 
which combines lexicographic bounds with traditional set bounds. They demonstrated that 
a set solver based upon lexicographic bounds consistency techniques produced stronger 
propagation than a traditional set bounds solver, although this came at a substantial cost 
in propagation performance. However, given that the use of the ROBDD representation 
leads to a performance improvement in the case of set bounds propagation, it is worth 
investigating the performance of an ROBDD-based lexicographic bounds propagator. 

The lexicographic bounds of a domain can be very compactly represented as an ROBDD. 
Like set bounds, an ROBDD representing the upper or lower lexicographic bounds of a 
domain is an ROBDD of size O(N), and so is the combination. Because these ROBDDs are 
compact this hopefully leads to fast propagation. Moreover, given an ROBDD domain, it is 
very easy to extract the lexicographic bounds on that domain in a single pass. By a process 
analogous to the construction of the bounds and cardinality propagator, we can use a split 
domain propagator combined with a function which determines the lexicographic bounds 
of a domain to construct a highly efficient lexicographic bounds solver. 

We define two functions bdd_lex_lower_bound and bdd_lex_upper_bound. These functions, 
given an ROBDD representing the domain D of a variable V together with a list of Boolean 
variables B corresponding to the bits of V, return the lower and upper lexicographic bounds 
(respectively) of D. bdd_lex_lower_bound is implemented as shown in Figure 15 (the imple- 
mentation of bdd_lex_upper_bound is very similar). 

We define: 

bddJex_bounds(D, B) = bddJexJower_bound(D, B) A bdd Jex_upper_bound(D, B) 

The split set and lexicographic bounds propagator can then be implemented exactly as 
in Equation (7), using bdd_lex_bounds in place of bdd_card_bounds. 

Example 6.5. Consider lexicographic bounds propagation on the constraint c = |s| = 2 
where s C {1,2,3,4,5}. The ROBDD B{c) is shown in Figure 3(b) and initial domain 
D(s) = [0, {1,2,3,4,5}]. Then 5 S = B(c) and (3 S = 1 so the call to bdd Jex_bounds cal- 
culates the lower bounds lexicographic ROBDD shown in Figure 16(a), the upper bounds 
lexicographic ROBDD shown in Figure 16(b), and final answer the conjunction shown in 
Figure 16(c). Note that we have lost some information relative to the original cardinality 
ROBDD. 

As observed in Section 5.3, the lexicographic ordering for set variables actually corre- 
sponds to a numeric ordering on integer variables, so a pure lexicographic bounds propagator 
would also be coincidentally an integer bounds propagator. 

Experimental results for the lexicographic bounds propagator are given in Section 7. 
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bddJex_lower_bound(Z), B): 

if (\B\ = or D = 1) return 1 

if (D = 0) error 

n(v, t, e) := D 

(b 1 ,...,b n ) := B 

if (6i y v) error 

else if (b\ = v and e = 0) 

return b\ A bdd_lexJower_bound(t, (62, • • • , b n )) 
else 

if (61 = v) r := bdd_lex_lower_bound(e, (62, . . . , b n )) 
else r := bdd_lex_lower_bound(D, (62, • • • , b n )) 
return 61 V (-161 A r) 



Figure 15: An algorithm to extract the lower lexicographic bounds on the domain of a set 
variable represented by an ROBDD D, where B is the vector of (non-fixed) bits 
in the Boolean representation of the set variable 




7. Experiments 

We have implemented a set solver using the ideas described in this paper. Our implemen- 
tation is written primarily in Mercury (Somogyi, Henderson, & Conway, 1996), using an 
interface to the C language ROBDD library CUDD as a platform for ROBDD manipula- 
tions (Somenzi, 2004). The ROBDD library is effectively treated as a black box. 
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We used our solver to implement a series of standard constraint benchmarks as described 
below. Many of these problems are in the CSPLib library of constraint satisfaction problems 
(Gent, Walsh, & Selman, 2004). For the purposes of comparison, we also implemented 
our benchmarks using the ic.sets library of ECL*PS e v5.6 2 and the finite sets library of 
Mozart vl.3.0. We conducted all of our experiments on a cluster of 8 identical 2.8GHz 
Pentium 4 machines, each with 1 Gb of RAM and running Debian GNU/Linux 3.1. All 
three solvers were limited to 1 Gb of memory to minimise swapping. All experiments were 
repeated three times, and the lowest time out of the three runs was taken as the result. 

7.1 Steiner Systems 

A commonly used benchmark for set constraint solvers is the calculation of small Steiner 
systems. A Steiner system S(t, k, N) is a set X of cardinality N and a collection C of 
subsets of X of cardinality k (called "blocks" ) , such that any t elements of X are in exactly 
one block. Steiner systems have been extensively studied in combinatorial mathematics. 
If t = 2 and k = 3, we have the so-called Steiner triple systems, which are often used as 
benchmarks (Gervet, 1997; Azevedo, 2002; Miiller, 2001). Any Steiner system must have 
exactly m = (^)/( k t ) blocks (Theorem 19.2 of van Lint k Wilson, 2001). 

It is natural to model a Steiner system using set variables si,...,s m , where each set 
variable corresponds to a single block, subject to the following constraints: 

TO 

A(N = fc ) (8) 
1=1 

m—1 m 

A f\ f\ (Bitij tiij = Si (~1 sj A \uij\ <{t- 1)) A (si < sj) (9) 

i=l j=i+l 

The lexicographic ordering constraint Sj < Sj has been added to remove symmetries from 
the problem formed by permuting the blocks. 

A necessary condition for the existence of a Steiner system is that (^_T/) / {^Zf) is an 
integer for all i G {0, 1, ...,£ — 1} (van Lint & Wilson, 2001); we say a set of parameters 
(t, k, N) is admissible if it satisfies this condition. In order to choose test cases, we ran 
each solver on every admissible set of (t, k, N) values for N < 32. Results are shown for 
every test case that at least one solver was able to solve within a time limit of 10 minutes, 
"x" denotes abnormal termination due to exceeding an arbitrary limit on the maximum 
number of ROBDD variables imposed by the CUDD package, while " — " denotes failure to 
complete a testcase within the time limit. 

In all cases we use a sequential variable-ordering heuristic, and a value-ordering heuristic 
that chooses the largest unfixed value within a variable's domain. Once a variable and 
an unfixed value have been chosen for labeling, either that value is a member of the set 
represented by the variable or it is not. The order in which we try the two alternatives 
has a significant effect on the performance of the solver. In this case we elect to choose the 
element-not-in-set option first. 

2. Very recently a new sets library Cardinal has been added to ECL l PS e which supports better cardinality 
reasoning. Unfortunately we cannot directly compare with it on these benchmarks since it does not 
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Table 2: Results for Steiner systems, with split constraints and intermediate variables. 



Testcase 


ECL*PS e 


Mozart 


Bounds 


Domain 


LU+R 


LU+Lex 


LU+Card 




Time 


Fails 


Time 


Fails 


Time 


Fails 


Time 


Fails 


Time 


Fails 


Time 


Fails 


Time 


Fails 




/s 




/s 








/* 




/" 




/* 




/s 




S(2,3,7) 


0.3 


10 


0.1 


21 


<0.1 


10 


0.1 





0.1 





0.1 


4 


<0.1 


2 


S(3,4,8) 


0.5 


21 


0.1 


52 


0.1 


21 


0.4 





0.4 





0.4 


4 


0.1 


4 


S(2,3,9) 


7.7 


1394 


1.0 


5102 


1.6 


1394 


1.0 


100 


1.3 


100 


3.3 


421 


2.4 


1072 


S(2,4,13) 


1.8 


313 


0.4 


1685 


0.6 


313 


1.7 


32 


1.5 


32 


2.1 


127 


0.7 


157 


S(2,3,15) 


3.6 


65 


0.5 


354 


2.2 


65 


20.2 





19.6 





20.4 


127 


3.3 


41 


S(3,4,16) 


67.5 


289 






X 


X 


X 


X 


X 


X 


X 


X 


X 


X 


S(2,5,21) 


3.2 


421 


0.4 


668 


2.3 


421 


110.2 





59.8 





21.5 


139 


2.6 


124 


S(3,6,22) 


49.7 


1619 






X 


X 


X 


X 


X 


X 


X 


X 


X 


X 



Table 3: Results for Steiner systems, with merged constraints and no intermediate variables. 



Testcase 


Bo 


unds 


Domain 


LI 


r+R 


LU+Lex 


LU+Card 




Time 


Fails 


Time 


Fails 


Time 


Fails 


Time 


Fails 


Time 


Fails 




/s 




h 




/s 




h 




h 




S(2,3,7) 


<0.1 


8 


<0.1 





<0.1 





<0.1 





<0.1 





S(3,4,8) 


0.1 


18 


0.1 





0.1 





0.1 





0.1 





S(2,3,9) 


0.2 


325 


0.1 


9 


0.1 


9 


0.1 


11 


0.2 


113 


S(2,3,13) 






109.2 


24723 


144.6 


24723 


518.1 


30338 






S(2,4,13) 


0.1 


157 


0.1 





0.1 





0.4 


11 


0.1 


27 


S(2,3,15) 


0.4 


56 


1.3 





1.4 





2.9 





0.7 


32 


S(2,4,16) 


421.4 


522706 


0.6 


15 


0.6 


15 


2.5 


16 


577.0 


209799 


S(2,6,16) 






80.7 


15205 


82.7 


15205 










S(3,4,16) 


9.7 


274 


548.7 





485.3 





428.9 





18.4 


162 


S(2,5,21) 


0.5 


413 


1.4 





1.4 





32.9 





0.6 


116 


S(3,6,22) 


8.3 


1608 














12.7 


381 


S(2,3,31) 


23.3 


280 














48.6 


224 



Table 4: All-solutions results on Steiner systems. " — " denotes failure to complete a test 
case within one hour 







ECL*PS e 


Bounds 


Domain 


LU+R 


LU+Lex 


LU+Card 


Problem 


Solns. 


time 

/* 


fails 


time 

/s 


fails 


time 

1* 


fails 


time 

/s 


fails 


time 

1* 


fails 


time 

/s 


fails 


S(2,3,7) 


30 


16.3 


3,015 


0.2 


537 


0.1 


47 


0.1 


47 


0.2 


76 


0.3 


267 


S(3,4,8) 


30 






726.5 


610271 


1.4 


492 


2.1 


492 


22.4 


3248 


951.2 


431801 


S(2,3,9) 


840 






398.1 


391691 


23.4 


16794 


37.8 


16794 


110.4 


29133 


593.3 


224131 


S(2,6,16) 













80.9 


15205 


83.0 


15205 
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In order to compare the raw performance of the various solvers, irrespective of any 
modeling advantages of the ROBDD-based solver, we performed experiments using a model 
of the problem which contains only primitive constraints and makes use of intermediate 
variables. This "split" model contains m unary constraints corresponding to Equation (8) 
and 3m(m— 1)/2 binary constraints corresponding to Equation (9) (containing intermediate 
variables Uij). The same model was used in the ECI/PS 6 , Mozart and ROBDD-based 
solvers, permitting a direct comparison of propagation performance. Results for this model 
are shown in Table 2. In particular, observe that with this model ECI/PS 6 and the ROBDD- 
based bounds solver produce the same number of failures, demonstrating that the search 
spaces explored by the two solvers are identical. 

One of the main strengths of the ROBDD-based modeling approach is that it gives us 
the freedom to merge arbitrary constraints and existentially quantify away intermediate 
variables, allowing us to model set constraint problems more efficiently In the case of 
Steiner systems, this allows us to model the problem as m unary constraints corresponding 
to Equation (8), and just m(m — l)/2 binary constraints ipij of the form ipij = (|sj fl Sj\ < 
(t — 1)) A (si < Sj). The binary constraints do not contain the intermediate variables 
Uij — they are not required since they can be existentially quantified out of the ROBDD 
representation. Experimental results for this revised model are shown in Table 3. 

In all cases the revised model propagates much more strongly than the original model, 
leading to a substantial decrease in solution time. In addition, the decrease in the number 
of set variables required permits the solution of larger test cases. Clearly it is beneficial to 
remove intermediate variables and merge constraints. 

Despite weaker propagation the ROBDD bounds solver is often the fastest method of 
finding a single solution to a Steiner System. In order to determine whether this was due 
to the efficiency of the solver, or whether the solver was just "lucky" in finding a solution 
quickly, we also ran experiments to find all solutions of the Steiner systems up to reordering 
of the blocks. The results for all test cases that at least one of the solvers was able to solve 
within a time limit of one hour are shown in Table 4. In all cases the reduction both of 
time and number of fails demonstrate the superiority of the propagation approaches based 
on domain consistency. 

7.2 Social Golfers 

Another problem often used as a set benchmark is the "Social Golfers" problem (problem 
probOlO of CSPLib). The aim of this problem is to arrange N = g x s golfers into g groups 
of s players for each of w weeks, such that no two players play together more than once. 
We can model this problem as a set constraint problem using a w x g matrix of set variables 
Vij, where 1 < i < w is the week index and 1 < j '• < g is the group index. 



support lexicographic ordering constraints. Testing without lexicographic orderings showed that it was 
about 2-3 times slower than LU + Card. 
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We use the following model of the problem: 

w w g 

/\ (partition^!, . . .,v ig )) A f\ f\ |%-| = s 

i=l i=l j=l 

w—1 w 

A A A \ Vik nv ji\^ iA /\ A m - v i i 

i,je{l,...,w} k,le{l,...,g} i=l j=i+l 

The partition < global constraint is a combined partitioning and lexicographical order- 
ing constraint, formed by merging the partition constraint of Section 4.2 with constraints 
imposing a lexicographic order on the variables. This constraint is trivial to construct using 
ROBDDs, but is not available in either ECL l PS e or Mozart. 3 

Results for each of the solvers are shown in Table 5 and Table 6. In the former table, 
a sequential smallest-element-in-set labeling strategy was used to enable a fair comparison 
of propagation performance, whereas in the latter table a first-fail labeling strategy was 
used in order to give a measure of the peak performance of each solver. For every test case 
in both tables, with a single exception (5-8-3), at least one of the ROBDD-based solvers 
performs equal or better to both ECL l PS e and Mozart. It should also be observed that 
when using a first-fail labeling strategy, the domain and split domain solvers are the only 
solvers able to solve every test case. 

There are several other features of the results that are worth noting. As in the case 
of Steiner systems, the ROBDD-based set bounds solver is often the fastest, despite weak 
propagation. Amongst the solvers with stronger propagation, the split domain solver is 
almost always faster than the original domain solver due to smaller domain sizes. It is, 
however, slower than the original domain solver in the presence of backtracking (due to 
the requirement to trail more values — in particular the propagator ROBDDs). The lexi- 
cographic bounds solver is almost as effective as the domain solvers in restricting search 
space, although it is usually outperformed by the domain and bounds solvers. 

7.3 Weighted Hamming Codes 

The problem of finding maximal Hamming Codes can be modelled as a set constraint 
problem. 

We define an ^-bit codeword to be a bit-string (or vector of Boolean values) of length /. 
Given two Z-bit codewords A and B, we define the Hamming distance d(A,B) between A 
and B to be the number of positions at which the two bit-strings differ. An (/, d)-Hamming 
Code is a set of Z-bit codewords such that the Hamming distance between any two codewords 
in the set is at least d. 

Given a codeword length I and the minimum Hamming distance d, the problem is to 
construct a Hamming code with the largest possible number of codewords. A variant of this 
problem, used as a benchmark by Sadler and Gervet (2004), has the additional requirement 
that each codeword have a fixed weight w, where the weight of a codeword is defined to be 

3. It would, of course, be possible to implement such a constraint in both ECL l PS e and Mozart, but 
such an implementation would be a fairly laborious process. A strength of the ROBDD-based modeling 
approach is that we can construct global constraints with no extra code. 
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Table 5: First-solution performance results on the Social Golfers problem, using a sequential 
"smallest-element-in-set" labeling strategy. Time and number of failures are given 
for all solvers. " — " denotes failure to complete a test case within 10 minutes. The 
cases 5-4-3, 6-4-3, and 7-5-5 have no solutions 



Problem 


ECL*PS e 


Mozart 


Bounds 


Domain 




LU+Lex 


LU+Card 


time 


fails 


time 


fails 


time 


fails 


time 


fails 


time 


time 


fails 


time 


fails 


■7/1— fl— Q 


/s 




1* 




h 




/s 




/s 






/s 




9^4 


7.6 


10468 


1.0 


7638 


0.1 


30 


0.1 







0.1 


3 


0.1 


5 


9 4 


49.2 


64308 


6.4 


42346 


0.6 


2036 


0.2 





n i 


0.7 


194 


0.4 


326 


9 7 4 


95.1 


114818 


10.9 


66637 


1.7 


4447 


0.4 





n 4 


2.3 


692 


1.6 


1608 
















2.0 





1 R 










^ ^ 4 


12.5 


14092 


2.5 


10311 


0.1 


30 


0.3 





u.o 


0.4 


3 


0.2 


5 


*\ (\ A 


76.3 


83815 


14.0 


51134 


1.6 


2039 


1.6 


o 


1 A 


2.2 


194 


1.9 


328 


3-7-4 


146.8 


146419 


27.3 


88394 


4.6 


4492 


8.9 





8.4 


7.6 


695 


5.1 


1629 


4-5-4 


14.1 


14369 


3.9 


10715 


0.2 


30 


0.8 





0.6 


0.7 


3 


0.4 


5 


4-6-5 










21.9 


12747 


118.6 





80.7 


19.3 


499 


32.0 


2122 


4-7-4 


169.3 


149767 


46.4 


90712 


8.7 


4498 






481.6 


14.1 


696 


10.5 


1632 


4-9-4 


27.3 


19065 


7.8 


12489 


2.6 


71 








22.2 


8 


8.9 


33 


5-4-3 










113.9 


63642 


28.6 


5165 


32.1 


88.9 


10210 


202.7 


50542 


5-5-4 


350.6 


199632 


217.7 


416889 


7.0 


2686 


3.8 


41 


2.3 


12 


313 


20.8 


1584 


5-7-4 










14.6 


4583 








23.7 


700 


17.5 


1683 


5-8-3 


5.0 


2229 


0.9 


1820 


1.1 


14 


9.2 





7.8 


3.2 


3 


2.1 


4 


6-4-3 










158.2 


61770 


20.3 


2132 


23.0 


60.8 


4506 


293.5 


49966 


6-5-3 


458.9 


240296 


287.4 


471485 


4.1 


1455 


1.8 


82 


1.5 


4.1 


202 


8.3 


1078 


6-6-3 


3.3 


1462 


1.0 


1462 


0.5 


5 


1.8 





1.4 


1.2 





0.7 


1 


7-5-5 














0.5 





0.4 


2.4 





1.3 


22 



the number of "1" bits that codeword contains. We will denote an instance of this problem 
by H(l,d,w). 

As proposed by Miiller and Miiller (1997), we can model this problem for n codewords 
using n set variables Si, where 1 < i < n . A codeword Cj corresponds to the characteristic 
function of the set Si, i.e. bit j is set in codeword Cj if and only if j G Si. The Hamming 
distance d(Ci, Cj) between two codewords Ci and Cj can be calculated from the associated 
sets Si and Sj thus: 

d{d, Cj) = i-\SinSj\- \{i, ...,i}\(SiU Sj)\ 



We can remove symmetries created by permuting the codewords by introducing lexico- 
graphic ordering constraints Si < Sj, for all 1 < i < j < n. The complete model of the 
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Table 6: First-solution performance results on the Social Golfers problem, using a first-fail 
"smallest-element-in-set" labeling strategy. Time and number of failures are given 
for all solvers. The cases 5-4-3, 6-4-3, and 7-5-5 have no solutions 



Problem 


ECL l PS e 


Mozart 


Bounds 


Domain 


T T T i T"» 

LU+R 


LU+Lex 


LU+Card 


time 


fails 


time 


fails 


time 


fails 


time 


fails 


time 


time 


fails 


time 


fails 


■7/1— fl— Q 






1* 








1* 




la 


h 




h 




9 ^ A 

Z.-')--t 


7.9 


10468 


1.1 


7638 


0.1 


30 


0.1 





yj, _l 


0.1 


3 


0.4 


447 


9 A 


51.3 


64308 


6.5 


42346 


0.6 


2036 


0.2 





n i 

u. ± 


0.7 


186 


3.8 


3820 


9 7 A 


99.9 


114818 


11.1 


66637 


1.7 


4447 


0.4 





n a 


2.0 


390 


4.6 


6424 


9 q c; 














1.9 















^ ^ A 


14.5 


14092 


2.6 


10311 


0.1 


44 


0.3 





u.o 


0.4 


5 


1.3 


481 


^ f\ A 


91.8 


83815 


15.1 


51134 


1.9 


2361 


1.1 







2.4 


209 


13.7 


3722 


1 7 A 


183.0 


146419 


28.7 


88394 


5.2 


5140 


1.9 





1 7 


6.2 


512 


15.6 


6067 


A ^ A 
-±- <)--t 


18.0 


14369 


4.1 


10715 


0.3 


47 


0.8 


o 


u.u 


0.7 


7 


2.6 


494 


4-6-5 










38.5 


19376 


62.5 





40.1 


24.6 


607 






4-7-4 


243.9 


149767 


49.6 


90712 


9.9 


5149 


6.5 





5.1 


10.6 


405 


29.5 


5717 


4-9-4 


40.9 


19065 


8.3 


12489 


2.7 


143 


152.0 





107.4 


22.6 


13 


12.5 


516 


5-4-3 










187.5 


103972 


23.2 


3812 


26.0 


91.4 


10422 


309.1 


72669 


5-5-4 


394.9 


199632 


224.5 


416889 


4.5 


2388 


2.5 


18 


1.7 


23.4 


776 


41.4 


4730 


5-7-4 










17.6 


5494 


18.2 





12.8 


16.9 


447 


47.1 


5473 


5-8-3 


5.4 


2229 


0.9 


1820 


1.1 


19 


4.5 





3.9 


3.4 


2 


2.7 


83 


6-4-3 










234.2 


90428 


14.6 


1504 


15.2 


54.0 


4013 


349.9 


59805 


6-5-3 


501.0 


240296 


294.2 


471485 


1.6 


495 


1.2 


34 


1.0 


58.2 


3787 


473.6 


68673 


6-6-3 


3.6 


1462 


1.0 


1462 






1.6 


7 


1.3 


5.1 


292 


1.0 


8 


7-5-3 














16.9 


528 


13.1 


288.8 


9829 






7-5-5 














0.5 





0.4 


2.4 





1.3 


22 



problem is: 

n 

f\\Si\=w (11) 

i=\ 

n—l n 

A A A (\S l nS,\ + \(S~US])\<l-d)A{S t <S j ) (12) 

i=l .7=1+1 

The constraints described by Equation (12) are implemented using a single ROBDD for 
each pair of i and j values. This is possible since we can model the integer addition and 
comparison operations as ROBDDs as described in Section 5, using the representation of 
the cardinality of a set variable as an integer expression as described in Section 5.5. 

In order to find an optimal solution, we initially set n to 1, and repeatedly solve instances 
of the problem, progressively incrementing n to find larger and larger codes. We prove 
optimality of a solution for n = k by failing to solve the problem for n = k + 1; we then 
know that the optimal value for n was k. 
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Table 7: Statistics for the 51 Weighted Hamming Code testcases that were solved by all of 
the solvers. 





Bounds 


Domain 


LU+R 


LU+Lex 


LU+Card 




time 


fails 


time 


fails 


time 


fails 


time 


fails 


time 


fails 




/s 




/s 








1* 




/s 




Mean 


17.7 


110034 


0.2 


210.7 


0.3 


210.7 


AA 


1886 


3.6 


6604 


Total 


903.3 




11.4 




14.2 




222.6 




184.6 




Minimum 


0.03 





0.03 





0.03 





0.03 





0.03 





25th Percentile 


0.03 


19 


0.03 





0.03 





0.03 


3 


0.03 





Median 


0.05 


196 


0.04 


2 


0.04 


2 


0.04 


20 


0.04 


2 


75th Percentile 


0.67 


5604 


0.06 


25 


0.06 


25 


0.18 


143.5 


0.09 


112 


Maximum 


415.55 


3021057 


4.16 


3740 


4.69 


3740 


113.9 


45667 


92.09 


211677 



Table 8: Weighted Hamming Code testcases that were solved by at least one but not all 
solvers 





Bounds 


Domain 


LU+R 


LU 


-l-Lex 


LU+Card 




time 


fails 


time 


fails 


time 


fails 


time 


fails 


time 


fails 




h 




h 




/* 




h 




/s 




H(8,4,4) 






1.6 


224 


1.9 


224 


8.5 


986 






H(9,4,3) 






11.3 


5615 


20.6 


5615 










H(9,4,6) 






25.4 


16554 


45.7 


16554 


321.9 


56599 






H(10,6,5) 






26.7 


16635 


29.6 


16635 


528.8 


169457 


428.3 


762775 



In order to assist with timing comparisons, we use the same set of problem instances 
as Sadler and Gervet (2004). We consider sets of values #(/, d, w) where I G {6, 7, 8, 9, 10}, 
d € {4, 6, 8, 10, 12}, and w € {3, 4, 5, 6, 7, 8}, with d < I and w < I (trivially there is at most 
one solution if d > I and none if w > I). There are 62 such testcases; some of these are almost 
identical — in particular testcases #(/, d, w) and H(l, d,l — w) have solutions that differ only 
through complementation of the bits. The ROBDD-based solver can solve all but seven test 
cases (namely #(9, 4, 4), #(9, 4, 5), #(10, 4, 3), #(10, 4, 4), #(10, 4, 5), #(10, 4, 6), #(10, 4, 7)), 
which in reality contains three pairs of mirror image testcases. 

Since there are too many results to list each testcase individually, performance statistics 
for the testcases that all the solvers were able to solve are shown in Table 7. Those cases 
that were solved by at least one but not all of the ROBDD-based solvers are shown in 
Table 8. 

Sadler and Gervet (2004) report results for this problem using set bounds and lexico- 
graphic bounds solvers implemented in ECL l PS e . Their solvers were able to solve 50 of the 
testcases with a time limit of 240 seconds for each testcase. Some individual testcases took 
upwards of 100 seconds. By contrast, the ROBDD-based domain solver is capable of solv- 
ing 55 testcases in 76.4 seconds in total. Clearly in this case enforcing domain consistency 
brings a considerable reduction in search space, leading to a highly efficient solver. 
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Moreover, the set bounds and lexicographic bounds solvers implemented using the 
ROBDD platform appear to perform considerably better than those of Sadler and Gervet. 
Due to the graphical nature of their presentation it is difficult to quantify this perfor- 
mance difference; exact figures are presented for only two testcases — namely H(9, 4, 7) and 
H(10, 6, 7). For the testcase H(9,4,7) the ROBDD set bounds and lexicographic bounds 
solvers each found and proved an optimal solution in 0.6 and 0.3 seconds respectively, 
compared with the > 240 and 167.1 seconds respectively quoted for the bounds and lexico- 
graphic solvers of Sadler and Gervet. Similarly for if (10, 6, 7) the ROBDD set bounds and 
lexicographic bounds solvers found and proved an optimal solution in 1.9 and 1.1 seconds 
respectively, as compared with > 240 and 98.5 seconds. Given this dramatic performance 
difference, it appears the additional modeling flexibility of the ROBDD-based solver pro- 
vides a substantial performance gain. 

It should be noted that the performance results reported by Sadler and Gervet were run 
on a slightly slower machine (a 2 Ghz Pentium 4 machine). Nonetheless, the contribution 
to the performance difference due to machine speed is dwarfed by the performance gap 
between the two solvers. 

7.4 Balanced Academic Curricula 

The Balanced Academic Curriculum problem (problem prob030 of CSPLib) involves plan- 
ning an academic curriculum over a sequence of academic periods in order to provide a 
balanced load in each period. 

A curriculum consists of m courses (1 < i < m) and n academic periods. Each course i 
has a set of prerequisites and an associated academic load U. Every course must be assigned 
to exactly one period. In any given period the total number of courses must be at least 
some minimum number c and at most a maximum number d. In addition, within any given 
period the total academic load of the courses must be at least some minimum load a and 
at most a maximum load b. Let R be the set of prequisite pairs where i and j are 

courses. Prerequisite relationships must be observed, so for each pair of courses £ R, 
if course i is scheduled in period p, then course j must be scheduled in a period strictly 
prior to p. 

A model of this problem using set variables and set constraints was proposed by Hnich 
et al. (2002), although no experimental results were presented for that model. In this 
"primal" model, we use one set variable Si per academic period. Each Si represents the 
set of courses assigned to academic period i, so k G Si if and only if course k is assigned to 
period i. We can model the problem using the following constraints: 

• (SI) Every course is taken exactly once: partition(Si, . . . , S n ) 

• (S2) The number of courses in a period is between c and d: 

yf =1 ((c<\Si\)A(\Si\<d)) 

• (S3) The total academic load in a period is between a and b: 

V" =1 ((a < wsum(S;, (h, . . . ,t m })) A (wsum(S , j, (h, . . .,t m )) < b)) 

• (S4) Prerequisites are respected: 

Vf =1 Vj =1 (A (c , p)ei? (pG^)-(c^^)) 
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In general the partition and wsum constraints of the primal set model can be very 
large, making domain propagation impractical. 

In addition to presenting results for the primal model, Choi et al. (2003) demonstrated 
that a substantial performance improvement can be obtained for this problem through the 
use of redundant models together with channelling constraints. 

We can obtain better results from a "dual" set model, where additional set variables 
are introduced to model each course. We define set variables Xi (1 < i < m) representing 
each course, such that if k G X{ then course i is assigned to period k. We can then define 
the following constraints: 

• (CX) Channelling constraints: V^V^i G Sj) <-» (j G X t ) 

• (XI) Each course may be assigned to at most one period: V^ 1 |Xj| = 1 

• (X2) Prerequisites are respected: V(i, j) G R lexlt(Xj, Xi) 

We define the dual model to consist of the constraints {S2, S3, CX, XI, X2}. Constraints 
SI and S4 are no longer required, since they are propagation redundant. 

Unfortunately, while the dual model performs better than the primal set model, it is still 
incapable of proving the optimality of the solutions to any of the problems. We can obtain 
stronger propagation by introducing redundant global constraints on the total academic 
load and number of courses of all academic periods, as suggested by Choi et al. (2003). 

For each academic period j define two integer variables lj, representing the academic 
load in period j and qj , representing the number of courses in period j . We can then define 
the following constraints on these new integer variables: 

• (CI1) Channeling constraints on If. V™ =1 wsum(Si, (t\, . . . , t m )) = U 

• (CI2) Channeling constraints on qf. V? =1 \Si\ = qi 

• (II) Range constraints on If. Vf =1 (a < U) A (k < b) 

• (12) Range constraints on qf. V" =1 (c < qi) A (% < d) 

• (13) All loads should be undertaken: Yli=i h = Yh=i ^ 

• (14) All courses must be taken: Ya=i % = m 

We then define the "hybrid dual" model consisting of the constraints {CX, XI, X2, 
CI1, CI2, II, 12, 13, 14}. Note that none of the original constraints S1-S4 remain. For 
completeness it is also possible to define a "hybrid primal" model consisting of constraints 
{SI, S4, CI1, CI2, II, 12, 13, 14}, although as before the large ROBDD sizes make domain 
propagation impractical. 

Experimental results for the Balanced Academic Curriculum problem are shown in Ta- 
ble 9. The timing results for the Hybrid Dual model are comparable to the best results 
obtained by Hnich et al. (2002) and Choi et al. (2003). The number of failures required to 
solve the problem is orders of magnitude smaller than the best results presented in either 
paper, emphasising the value of domain propagation in this case. These examples also show 
a case where the LU+Lex solver is competitive and clearly the most robust solver over the 
different models. 
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Table 9: Performance results for the Balanced Academic Curriculum Problem for the 8, 10 
and 12 period problems. The first column contains the number of periods, the 
second column contains the model type (HD=hybrid dual, D=dual, P=primal, 
HP=hybrid primal), and the third column contains the maximum load per period 
b. In all cases a sequential "smallest-element-in-set" labeling method was used 
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Steiner system 5(2, 3, 15) 




100 200 300 100 200 

Labelling step Labelling step 



Social Golfers problem 4-5-4 




Labelling step Labelling step 

Figure 17: Comparison of total domain and total ROBDD sizes with labeling step for set 
bounds, set domain, split set domain, set and lexicographic bounds, and set and 
cardinality bounds solvers on the Steiner System 5(2, 3, 15) and Social Golfers 
problem 4-5-4. Note that the y axes of the total domain size graphs have a 
logarithmic scale (base 2). 



7.5 Comparing the Propagation Performance of the Solvers 

It is instructive to compare the propagation performance of the various ROBDD-based 
solvers graphically. Here we pick two small examples, namely the Steiner System 5(2, 3, 15) 
and the Social Golfers problem 4-5-4, using the default labeling for each problem, and graph 
BDD and domain sizes against number of labeling steps and time. 
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Steiner system 5(2, 3, 15) Social Golfers problem 4-5-4 

Figure 18: Comparison of total domain size with time for set bounds, set domain, split set 
domain, set and lexicographic bounds, and set and cardinality bounds solvers 
on the Steiner System 5(2,3,15) and Social Golfers problem 4-5-4. Note that 
the y axes have a logarithmic scale (base 2). 



Figure 17 depicts how the total domain sizes and total ROBDD sizes (total number of 
internal nodes) vary with each labeling step for each of the solvers for the Steiner System 
5(2, 3, 15) and for the Social Golfers problem 4-5-4. In particular, observe that the domain 
and split-domain propagators are the most effective at reducing the domain size and thus 
restricting the search space, although maintaining domain consistency can be costly due to 
the size of the ROBDDs required for representing arbitrary domains. This can be seen in 
the ROBDD size against labeling step graph, most notably in the case of the Social Golfers 
problem 4-5-4. The lexicographic bounds solver is the next most effective in domain size 
reduction, and requires a smaller number of ROBDD nodes to store its bounds than the 
domain or cardinality bounds solvers. The bounds solver is clearly the weakest in terms of 
domain size reduction, and as opposed to the other solvers starts from a zero size domain 
representation and builds slowly to ROBDDs representing the answer. 

We are not only concerned with the strength of a propagator in the abstract, since 
the efficiency of the solver as a whole is dependent on both the propagation and labeling 
processes. Figure 18 depicts graphs of domain size against time for each of the propa- 
gators discussed in this paper. Here we see that the most effective propagator does not 
necessarily lead to the most efficient set solver. Despite comparatively weak propagation, 
the set and cardinality bounds solvers lead to the fastest reduction in domain size per unit 
time. Nonetheless, the experimental timing results given above demonstrate that for harder 
cases the additional cost of maintaining domain consistency or lexicographic bounds can be 
justified through the consequent reduction in search space. 
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Interestingly we can see from the graphs in Figure 18 that the initial domain reduction 
before labeling, which is the gap from time to where the line starts, can be a significant 
part of the computation. In particular, one of the weaknesses of the lexicographic bounds 
solver is the large time required to reach an initial fixpoint through lexicographic bounds 
propagation. For practical applications it might be preferable to implement a hybrid solver 
which utilises one of the other propagators to generate the initial domains, and uses the 
lexicographic bounds propagator during labeling. 

8. Conclusion 

We have demonstrated that ROBDDs form a highly flexible platform for constructing set 
constraint solvers. ROBDDs allow a compact representation of many set domains and 
set constraints, making them an effective basis for an efficient set solver. Since ROBDDs 
can represent arbitrary Boolean formulae we can easily conjoin and existentially quantify 
them, permitting the removal of intermediate variables and making the construction of 
global constraints trivial. We have demonstrated how to efficiently enforce various levels of 
consistency, including set domain, set bounds, cardinality bounds and lexicographic bounds 
consistency. Finally, we have presented experimental results that demonstrate that the 
ROBDD-based solver outperforms other common set solvers on a wide variety of standard 
set constraint satisfaction problems. 

No single set solver is uniformly better than the others. For many examples simple 
bounds propagation is the best approach, while in other cases, particularly when we ask 
for all solutions, domain consistency is preferable. There are also examples where lexico- 
graphic bounds or cardinality bounds are the best approach. The split-domain approach 
is somewhat disappointing, since it appears that in many cases the overhead of the more 
complicated calculation (Equation 6) is not rewarded in terms of smaller ROBDD sizes. 

In the future we plan to explore a robust general set constraint solver that dynamically 
chooses which level of consistency to maintain by examining how big the domain ROBDDs 
are becoming as search progresses. 
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